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Abstract 

We report on measurements of neutrino oscillation using data from the T2K long-baseline neu¬ 
trino experiment collected between 2010 and 2013. In an analysis of muon neutrino disappearance 
alone, we find the following estimates and 68% confidence intervals for the two possible mass 
hierarchies: 

Normal Hierarchy: sin 2 023 = 0.514 +qq5q and A = (2.51 ±0.10) x 10~ 3 eV 2 /c 4 
Inverted Hierarchy: sin 2 023 = 0.511 ± 0.055 and Amf 3 = (2.48 ± 0.10) x 10~ 3 eV 2 /c 4 

The analysis accounts for multi-nucleon mechanisms in neutrino interactions which were found to 
introduce negligible bias. 

We describe our first analyses that combine measurements of muon neutrino disappearance and 
electron neutrino appearance to estimate four oscillation parameters, |Am 2 |, sin 2 023, sin 2 0 i 3 , Scp , 
and the mass hierarchy. Frequentist and Bayesian intervals are presented for combinations of these 
parameters, with and without including recent reactor measurements. At 90% confidence level and 
including reactor measurements, we exclude the region 5cp = [0.15, 0.83]7t for normal hierarchy and 
$CP = [—0.08,1.09]7T for inverted hierarchy. The T2K and reactor data weakly favor the normal 
hierarchy with a Bayes Factor of 2.2. The most probable values and 68% ID credible intervals for 
the other oscillation parameters, when reactor data are included, are: 

sin 2 023 = 0.528^Q3g and | Am§ 2 | = (2.51 ± 0.11) x 10" 3 eV 2 /c 4 . 
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I. INTRODUCTION 


Neutrino oscillation was firmly established in the late 1990’s with the observation by the 
Super-Kamiokande (SK) experiment that muon neutrinos produced by cosmic ray interac¬ 
tions in our atmosphere changed their flavor [T]. Measurements from the Sudbury Neutrino 
Observatory a few years later, in combination with SK data, revealed that neutrino oscilla¬ 
tion was responsible for the apparent deficit of electron neutrinos produced in the Sun [2], In 
the most recent major advance, the T2K experiment M and reactor experiments |5H8] have 
established that all three neutrino mass states are mixtures of all three flavor states, which 
allows the possibility of CP violation in neutrino oscillation. This paper describes our most 
recent measurements of neutrino oscillation including our first results from analyses that 
combine measurements of muon neutrino disappearance and electron neutrino appearance. 

The Tokai to Kamioka (T2K) experiment [9] was made possible by the construction of 
the J-PARC high-intensity proton accelerator at a site that is an appropriate distance from 
the SK detector for precision measurements of neutrino oscillation. Protons, extracted from 
the J-PARC main ring, strike a target to produce secondary hadrons, which are focused 
and subsequently decay in-flight to produce an intense neutrino beam, consisting mostly of 
muon neutrinos. The neutrino beam axis is directed 2.5 degrees away from the SK detector, 
in order to produce a narrow-band 600 MeV flux at the detector, the energy that maximizes 
muon neutrino oscillation at the 295 km baseline. Detectors located 280 m downstream of 
the production target measure the properties of the neutrino beam, both on-axis (INGRID 
detector) and off-axis in the direction of SK (ND280 detector). 

T2K began operation in 2010 and was interrupted for one year by the Great East Japan 
Earthquake in 2011. The results reported in this paper use data collected through 2013, as 
summarized in Tab. [TJ With these data, almost 10% of the total proposed for the experiment, 
T2K enters the era of precision neutrino oscillation measurements. In 2014, we began to 
collect our Erst data in which the current in the magnetic focusing horns is reversed, so 
as to produce a beam primarily of muon anti-neutrinos. Future publications will report on 
measurements using that beam configuration. 

We begin this paper by describing the neutrino beamline and how we model neutrino 
production and interactions. We then summarize the near detectors and explain how we use 
their data to improve model predictions of neutrino interactions at the far detector. This 
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is followed by an overview of the far detector, how neutrino candidate events are selected, 
and how we model the detector response. Next, we describe the neutrino oscillation model, 
list the external inputs for the oscillation parameters, summarize the approaches used in the 
oscillation analyses, and characterize our main sources of systematic uncertainty. The final 
sections give detailed descriptions and results for the analysis of disappearance alone [ID] 
and for the joint analyses of disappearance and v e appearance. 

TABLE I: T2K data-taking periods and the protons on target (POT) used in the analyses 
presented in this paper. The maximum stable proton beam power achieved was 230 kW. 


Run Period Dates POT 

Run 1 Jan. 2010-Jun. 2010 0.32 x 10 20 

Run 2 Nov. 2010-Mar. 2011 1.11 x 10 20 

Run 3 Mar. 2012-Jun. 2012 1.58 x 10 20 

Run 4 Oct. 2012-May 2013 3.56 x 10 20 

Total Jan. 2010-May 2013 6.57 x 10 20 








II. NEUTRINO BEAMLINE 


The T2K primary beamline transports and focuses the 30 GeV proton beam extracted 
from the J-PARC Main Ring onto a 91.4-cm long graphite target. The secondary beamline 
consists of the target station, decay volume, and beam dump. The apparatus has been 
described in detail elsewhere [9] . 

The upstream end of the target station contains a collimator to protect the three down¬ 
stream focusing horns. The graphite target sits inside the first horn, and pions and other 
particles exiting the target are focused by these magnetic horns and are allowed to decay in 
the 96-m-long decay volume. Following the decay volume, protons and other particles that 
have not decayed are stopped in a beam dump consisting of 3.2 m of graphite and 2.4 m of 
iron, while muons above 5 GeV pass through and are detected in a Muon Monitor, designed 
to monitor the beam stability. With further absorption by earth, a beam of only neutrinos 
(primarily z/ M ) continues to the near and far detectors. 


A. Neutrino flux simulation 

The secondary beamline is simulated in order to estimate the nominal neutrino flux (in 
absence of neutrino oscillations) at the near and far detectors and the covariance arising 
from uncertainties in hadron production and the beamline configuration m We use the 
FLUKA 2008 package [I2l 113] to model the interactions of the primary beam protons and 
the subsequently-produced pions and kaons in the graphite target. As described below, we 
tune this simulation using external hadron production data. Particles exiting the target are 
tracked through the magnetic horns and decay volume in a GEANT3 [13] simulation using 
the GCALOR pA] package to model the subsequent hadron decays. 

In order to precisely predict the neutrino flux, each beam pulse is measured in the primary 
neutrino beamline. The suite of proton beam monitors consists of five current transformers 
which measure the proton beam intensity, 21 electrostatic monitors which measure the pro¬ 
ton beam position, and 19 segmented secondary emission monitors and an optical transition 
radiation monitor [[16] which measure the proton beam profile. The proton beam proper¬ 
ties have been stable throughout T2K operation, and their values and uncertainties for the 
most recent T2K run period, Run 4, are given in Tab. [TT] The values for other run periods 
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TABLE II: Summary of the estimated proton beam properties and their systematic errors 
at the collimator for the T2K Run 4 period. Shown are the mean position (X, Y), angle 
(X',Y'), width (cr), omittance (e), and Twiss parameter (cc) [17 ], 


X Profile Y Profile 
Parameter Mean Error Mean Error 


X, Y (mm) 

0.03 

0.34 

-0.87 

0.58 

X',Y’ (mrad) 

0.04 

0.07 

0.18 

0.28 

a (mm) 

3.76 

0.13 

4.15 

0.15 

e ( 7 r mm mrad) 

5.00 

0.49 

6.14 

2.88 

a 

0.15 

0.10 

0.19 

0.35 


have been published previously m The neutrino beam position and width stability is also 


monitored by the INGRID detector, and the results are given in Sec. IV A 


To improve the modeling of hadron interactions inside and outside the target, we use 
data from the NA61/SHINE experiment [T8j [T9] collected at 31 GeV/c and several other 
experiments [2UR22] . The hadron production data used for the oscillation analyses described 
here are equivalent to those used in our previous publications mm, including the statistics- 
limited NA61/SHINE dataset taken in 2007 on a thin carbon target. The NA61/SHINE 
data analyses of the 2009 thin-target and T2K-replica-target data are ongoing, and these 
additional data will be used in future T2K analyses. We incorporate the external hadron 
production data by weighting each simulated hadron interaction according to the measured 
multiplicities and particle production cross sections, using the true initial and final state 
hadron kinematics, as well as the material in which the interaction took place. The predicted 
flux at SK from the T2K beam is shown in Fig. [T| 


B. Neutrino flux uncertainties 

Uncertainty in the neutrino flux prediction arises from the hadron production model, 
proton beam profile, horn current, horn alignment, and other factors. For each source of 
uncertainty, we vary the underlying parameters to evaluate the effect on the flux prediction 


in bins of neutrino energy for each neutrino flavor Enj. Table |III| shows the breakdown for 


10 









TABLE III: Contributions to the systematic uncertainties for the unoscillated and v e 
flux prediction at SK, near the peak energy and without the use of near detector data. 
The values are shown for the (y e ) energy bin 0.6 GeV < E v < 0.7 GeV (0.5 GeV 

< E u < 0.7 GeV). 


Error source 

Uncertainty in 

SK flux near peak (%) 


u ii 

v e 

Beam current normalization 

2.6 

2.6 

Proton beam properties 

0.3 

0.2 

Off axis angle 

1.0 

0.2 

Horn current 

1.0 

0.1 

Horn field 

0.2 

0.8 

Horn misalignment 

0.4 

2.5 

Target misalignment 

0.0 

2.0 

MC statistics 

0.1 

0.5 

Hadron production 



Pion multiplicities 

5.5 

4.7 

Kaon multiplicities 

0.5 

3.2 

Secondary nucleon multiplicities 

6.9 

7.6 

Hadronic interaction lengths 

6.7 

6.9 

Total hadron production 

11.1 

11.7 

Total 

11.5 

12.4 


the and v e flux uncertainties for energy bins near the peak energy. 

The largest uncertainty from beam monitor calibrations arises in the beam current mea¬ 
surement using a current transformer, but its effect on the oscillation analyses is reduced 
through the use of near detector data. The remaining uncertainties due to the uncertain 
position and calibration of the other beam monitors are significantly smaller. As described 
in Sec. IV A the neutrino beam direction is determined with the INGRID detector, and 
therefore the assigned uncertainty on the off-axis angle comes directly from the INGRID 
beam profile measurement. To account for the horn current measurement that drifts over 
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time and a possible scale uncertainty, 5 kA is assigned as a conservative estimate of the 
horn current error. In the flux simulation, the horn magnetic field is assumed to have a 1 jr 
dependence. Deviations from this field, measured using a Hall probe, are used to define the 
uncertainty of the horn field. Horn and target alignment uncertainties come from survey 
measurements. 

Systematic uncertainties in modeling particle multiplicities from hadronic interactions 
come from several sources: experimental uncertainties in the external data, the uncertain 
scaling to different incident particle momenta and target materials, and extrapolation to 
regions of particle production phase space not covered by external data HU- The overall 
uncertainty is described by calculating the covariance of the pion, kaon, and secondary 
nucleon multiplicities and their interaction lengths. 

The systematic errors on the flux at SK, without applying near detector data, are 
shown in bins of neutrino energy in Fig. [2j The dominant source of uncertainty is from 
hadron production. 

For analyses of near and far detector data, the uncertainties arising from the beamline 
configuration and hadron production are propagated using a vector of systematic parameters, 
b, which scale the nominal flux in bins of neutrino energy, for each neutrino type (u e , z/ M , u e , 
Vfj) at each detector (ND280 and SK). The energy binning for each neutrino type is shown 
in Fig. [lj The covariance for these parameters is calculated separately for each T2K run 
period given in Tab. |TJ and the POT-weighted average is the flux covariance, 14, used by the 
near detector and oscillation analyses. We define b n and b s as the sub-vector elements of b 
for ND280 and SK. It is through the covariance between b n and b s that the near detector 
measurements of events constrain the expected unoscillated far detector and v e event 
rates in the oscillation analyses. 
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FIG. 1: The T2K unoscillated neutrino flux prediction at SK is shown with bands 
indicating the systematic uncertainty prior to applying near detector data. The flux in the 
range 8 GeV < E v < 30 GeV is simulated but not shown. The binning for the vector of 
systematic parameters, b , for each neutrino component is shown by the four scales. The 
same binning is used for the ND280 and SK flux systematic parameters, b n and b s . 
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FIG. 2: Fractional systematic error on the flux at SK arising from the beamline 
configuration and hadron production, prior to applying near detector data constraints. 
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III. NEUTRINO INTERACTION MODEL 


Precision neutrino oscillation measurements rely on having an accurate neutrino interac¬ 
tion model. The model is used to evaluate the selection efficiencies of the different signal 
and background interactions as well as the estimate of the neutrino energy from the detected 
final state particles. Finally, the model forms the basis to account for differences in the pre¬ 
dicted neutrino cross sections between different T2K detectors due to their different target 
nuclei compositions. All of these factors and their uncertainties are incorporated into the 


model for the T2K experiment through a set of systematic parameters x listed in Tab. VI] 
and their covariance V x . 

This section describes the interaction model in NEUT, the primary neutrino interaction 
generator used by T2K, explains how we use data from external experiments to provide 
initial constraints on the model before fitting to T2K data, discusses remaining uncertainties 
not constrained by external data sources, and discusses uncertainties based on differences 
between the NEUT model and those found in other interaction generators. 


A. Neutrino Interaction Model 

The interaction model used in this analysis is NEUT [ 23 ] version 5.1.4.2, which models 
neutrino interactions on various nuclear targets over a range of energies from ~100 MeV to 
~100TeV. NEUT simulates seven types of charged current (CC) and neutral current (NC) 
interactions: (quasi-)elastic scattering, single pion production, single photon production, 
single kaon production, single eta production, deep inelastic scattering (DIS), and coherent 
pion production. Interactions not modeled in this version of NEUT include, but are not 
limited to, multi-nucleon interactions in the nucleus [24, 25], and neutrino-electron scattering 
processes. 

The Llewellyn Smith model [25] is used as the basis to describe charged current quasi¬ 
elastic (CCQE) and neutral current elastic scattering (NCEL) interactions. In order to 
take into account the fact that the target nucleon is in a nucleus, the Relativistic Fermi 
Gas (RFG) model by Smith and Moniz [27, El] is used. The model uses dipole axial form 
factors and the vector form factors derived from electron scattering experiments [29]. The 
default quasi-elastic axial mass, is 1.21 GeV/c 2 and the default Fermi momenta for the 
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two dominant target nuclei carbon and oxygen are 217 MeV/c and 225 MeV/c, respectively. 
Appropriate Fermi momenta, pp, and binding energies, E B , are assigned to the other target 
nuclei. 

The Rein and Sehgal model [3U]J is used to simulate neutrino-induced single pion produc¬ 
tion. The model assumes the interaction is split, into two steps as follows: v + N —y l + N*, 
N* —y tt + N', where N and N' are nucleons, i is an outgoing neutrino or charged lepton, 
and N * is the resonance. For the initial cross section calculation, the amplitude of each 
resonance production is multiplied by the branching fraction of the resonance into a pion 
and nucleon. Interference between 18 resonances with masses below 2 GeV/c 2 are included 
in the calculation. To avoid double counting processes that produce a single pion through 
either resonance or DIS in calculating the total cross section, the invariant hadronic mass 
W is restricted to be less than 2 GeV/c 2 . The model assigns a 20% branching fraction for 
the additional delta decay channel that can occur in the nuclear medium, A + N —y N + IV, 
which we refer to as pion-less delta decay (PDD). Since the Rein and Sehgal model provides 
the amplitudes of the neutrino resonance production, we adjust the NEUT predictions for 
the cross sections of single photon, kaon, and eta production by changing the branching 
fractions of the various resonances. 

The coherent pion production model is described in [3T]. The interaction is described as 
u + A^i + n + X, where A is the target nucleus, i is the outgoing lepton, tt is the outgoing 
pion, and X is the remaining nucleus. The CC component of the model takes into account 
the lepton mass correction provided by the same authors [ 32] . 

The DIS cross section is calculated over the range of W > 1.3 GeV/c 2 . The structure 
functions are taken from the GRV98 parton distribution function [3.3] with corrections pro¬ 
posed by Bodek and Yang [M] to improve agreement with experiments in the low-Q 2 region. 
To avoid double counting single pion production with the resonance production described 
above, in the region W < 2 GeV/c 2 the model includes the probability to produce more than 
one pion only. For W > 2 GeV/c 2 , NEUT uses PYTHlA/JetSet [35] for hadronization while 
for W < 2 GeV/c 2 it uses its own model. 

Hadrons that are generated in a neutrino-nucleus interaction can interact with the nucleus 
and these final state interactions (FSI) can affect both the total number of particles observed 
in a detector and their kinematics. NEUT uses a cascade model for pions, kaons, etas, 
and nucleons. Though details are slightly different between hadrons, the basic procedure 
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is as follows. The starting point for the cascade model is the neutrino interaction point 
in the nucleus based on a Woods-Saxon density distribution [35] except in DIS, where a 
formation zone is taken into account. The hadron is moved a small distance and interaction 
probabilities for that step are calculated. The interaction types include charge exchange, 
inelastic scattering, particle production, and absorption. If an interaction has occurred, then 
the kinematics of the particle are changed as well as the particle type if needed. The process 
is repeated until all particles are either absorbed or escape the nucleus. 


B. Constraints From External Experiments 

To establish prior values and errors for neutrino-interaction systematic parameters x and 
constrain a subset for which ND280 observables are insensitive, neutrino-nucleus scattering 
data from external experiments are used. 

The datasets external to T2K come from two basic sources: pion-nucleus and neutrino- 
nucleus scattering experiments. To constrain pion-nucleus cross section parameters in the 
NEUT FSI model, pion-nucleus scattering data on a range of nuclear targets are used. The 
most important external source of neutrino data for our interaction model parameter con¬ 
straints is the MiniBooNE experiment [37]. The MiniBooNE flux [38] covers an energy range 
similar to that of T2K and as a 47 t detector like SK has a similar phase space acceptance, 
meaning NEUT is tested over a broader range of Q 2 than current ND280 analyses. 


1. Constraints From Pion-Nucleus Scattering Experiments 

To evaluate the uncertainty in the pion transport model in the nucleus, we consider 
the effects of varying the pion-nucleus interaction probabilities via six scale factors. These 
scale factors affect the following processes in the cascade model: absorption ( X FSABS ), low 
energy QE scattering including single charge exchange (x FS ® E ) and low energy single charge 
exchange (SCX) (x FSCX ) in a nucleus, high energy QE scattering (x FS ® EH ), high energy 
SCX (x FSCXH ), and pion production (x FSINEL ). The low (high) energy parameters are used 
for events with pion momenta below (above) 500 MeV/c with the high energy parameters 
explicitly given and the remaining parameters all low energy. The simulation used to perform 
this study is similar to the one in [39]. The model is fit to a range of energy-dependent cross 
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TABLE IV: NEUT FSI parameters, x FSI , that scale each interaction cross section. Shown 
are the best-fit and the maximum and minimum scaling values from the 16 parameter sets 

taken from the 6-dimensional lcr surface. 



X FSQE 

X FSQEH 

X FSINEL 

X FSABS 

X FSCX 

x FSCXH 

Best Fit 

1.0 

1.8 

1.0 

1.1 

1.0 

1.8 

Maximum 

1.6 

2.3 

1.5 

1.6 

1.6 

2.3 

Minimum 

0.6 

1.1 

0.5 

0.6 

0.4 

1.3 


sections comprising nuclear targets from carbon to lead 


The best-fit scale factors 


for these parameters are shown in Tab. IV as well as the maximum and minimum values for 
each parameter taken from 16 points on the lcr surface of the 6-dimensional parameter space. 
The parameter sets are used for assessing systematic uncertainty in secondary hadronic 
interactions in the near and far detectors, as discussed in Secs. [V)l and [VTjb, respectively. 


2. Constraints From MiniBooNE CCQE Measurements 


To constrain parameters related to the CCQE model and its overall normalization, we fit 
the 2D cross-section data from MiniBooNE [67], binned in the outgoing muon kinetic energy, 
T /; , and angle with respect to the neutrino beam direction, 9^. The NEUT interactions 
selected for the fit are all true CCQE interactions. Our fit procedure follows that described 
by Juszczak et al. [68], with the y 2 defined as 


V(MfU) = 



pEpMaCF 

Api 


+ 


A- 1 - 1 
AA 


( 1 ) 


where the index i runs over the bins of the (T /t , cos 9^) distribution, pf p is the measured (pre¬ 
dicted) differential cross section, A pi is its uncertainty, A is the CCQE normalization, and 
AA is the normalization uncertainty, set at 10.7% by MiniBooNE measurements. The main 
difference from the procedure in [6S] is that we include (X/,cos6Q) bins where a large per¬ 
centage of the events have 4-momentum transfers that are not allowed in the RFG model. 
We find M® E = 1-64 ± 0.03 GeV/c 2 and A = 0.88±0.02 with y% n /DOF = 26.9/135. It 
should be noted that MiniBooNE does not report correlations, and without this informa¬ 
tion assessing the goodness-of-ht is not possible. To take this into account, we assign the 
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uncertainty to be the difference between the fit result and nominal plus the uncertainty on 
the fit result. The M® E fit uncertainty is set to 0.45 GeV/c 2 , which covers (at 1 standard 
deviation) the point estimates from our fit to the MiniBooNE data, the K2K result [59] and 
a world deuterium average, 1.03 GeV/c 2 [TO]. The normalization uncertainty for neutrinos 
with E u < 1.5 GeV, xf E , is set to 11%, the MiniBooNE flux normalization uncertainty, 
since most of the neutrinos from MiniBooNE are created in this energy range. 


3. Constraints From MiniBooNE Inclusive n Measurements 

To constrain single pion production parameter errors, we use published MiniBooNE dif¬ 
ferential cross-section datasets for CC single 7r° production (CC 17 T 0 ) [71], CC single tt + 
production (CCl7r + ) [72], and NC single 7r° production (NCl7r°) [72]. Because the modes 
are described by a set of common parameters in NEUT, we perform a joint fit to all three 
data sets. 

The selection of NEUT simulated events follows the signal definition in each of the Mini¬ 
BooNE measurements. For the (CCl7r°, CCl7r + , NCl7r°) selections, the signals are defined 
as (z/ M , zz M , u) interactions with (1,1,0) /i~ and exactly one (7r°,7r + ,7r°) exiting the target 
nucleus, with no additional leptons or mesons exiting. In all cases, there is no constraint on 
the number of nucleons or photons exiting the nucleus. 

We consider a range of models by adjusting 9 parameters shown in Tab. 0 M RES 

is 

the axial vector mass for resonant interactions, which affects both the rate and Q 2 shape 
of interactions. The “W shape” parameter is an empirical parameter that we introduce 
in order to improve agreement with NCl7T° |p^o | data. The weighting function used is a 
Breit-Wigner function with a phase space term: 

r(w ' s) = a ' (w-w°r + sy 4 ' p(w ' m ” m "> (2) 

where S is the “W shape” parameter, W 0 = 1218 MeV/c, P(W\ m n , m^) is the phase space 
for a two-body decay of a particle with mass W into particles with masses m n and mjy, and a 
is a normalization factor calculated to leave the total nucleon-level cross section unchanged 
as S is varied. The nominal values of S and W$ come from averages of fits to two W 
distributions of NEUT interactions, one with a resonance decaying to a neutron and 7r + and 
the other with it decaying to a proton and 7T°. The “CCOther shape” parameter, x ccoth , 


19 



modifies the neutrino energy dependence of the cross section for a combination of CC modes, 
as described in Sec. Ill C[ along with the remaining parameters that are normalizations 
applied to the NEUT interaction modes. Simulated events modified by x CCOth constitute 
a small fraction of the selected samples. As a result, the data have minimal power to 
constrain this parameter and likewise for the NCl7r + , NC coherent pion, and NCOther 
normalization parameters, x NCln± , x NCcohn , and x NCOth , respectively. The T2K oscillation 
analyses are insensitive to these poorly determined parameters, and an arbitrary constraint 
is applied to stabilize the fits. In our external data analysis the NC coherent normalization 
cannot be constrained independently of the NCl7T° normalization, x NC1 ’’ r °, because there 
is no difference in the |p^o | spectrum between the two components. The errors given in 
Tab. |V| also include the variance observed when refitting using the 16 FSI ler parameter sets 
and scaling the errors when fitting multiple datasets following the approach of Maltoni and 
Schwetz m ■ The “W shape” nominal prior is kept at the default of 87.7 MeV/c 2 and in 
the absence of reported correlations from MiniBooNE, the uncertainty is estimated as the 
difference between the best fit and default values. The correlations between M EES , x ECln , 
and x NClna are given in Table 


VI 


C. Other NEUT Model Parameters 

The remaining uncertainties are in the modeling of the CC resonant, CCDIS, NC resonant 
charged pion, CC and NC coherent pion, anti-neutrino, as well as v e CCQE interactions. 
An additional set of energy-dependent normalization parameters is added for CCQE and 
CC17T interactions. Finally, a normalization parameter for the remaining NC interactions is 
included. 

The CCOther shape parameter, x CCOih \ accounts for model uncertainties for CCDIS and 
resonant interactions where the resonance decays to a nucleon and photon, kaon, or eta. 
The nominal interaction model for these interactions is not modified. From MINOS HE 
the uncertainty of their cross section measurement at 4 GeV, which is dominated by CCDIS, 
is approximately 10%. Using this as a reference point, the cross section is scaled by the 
factor (1 + x ccoth / E u ) where E u is the neutrino energy in GeV. The nominal value for 
xpeoth j g q an( j p as a i a constraint of 0.4. 

Normalization parameters are included for both CC and NC coherent pion interactions, 
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TABLE V: Parameters used in the single pion fits and their results from fitting the 
MiniBooNE data. Those with an arbitrary constraint applied have their la penalty term 
shown. M EES , x < p c ' l7r , and x NCln ° fit results and their covariance are used in subsequent 

analyses. 


units Nominal value Penalty Best fit Error 


M RES 

GeV/c 2 

1.21 


1.41 

0.22 

W shape 

MeV/c 2 

87.7 


42.4 

12 

gCCcoh-TT 


1 


1.423 

0.462 

™CCl7T 

X 1 


1 


1.15 

0.32 

x CCOth 


0 

0.4 

0.360 

0.386 

r^NCcohn 


1 

0.3 

0.994 

0.293 

x NClir° 


1 


0.963 

0.330 

X NC1t± 


1 

0.3 

0.965 

0.297 

x NCOth 


1 

0.3 

0.987 

0.297 


TABLE VI: Correlation between M EES , xf ci7r , and x NCl7T ° 



m r es 

™CCl7T 

X 1 

x NC1tt° 

M ees 

1 

-0.26 

-0.30 

„CCln 

X 1 

-0.26 

1 

0.74 

x NCln ° 

-0.30 

0.74 

1 


xpCcoh-K an( ] x NCcoh-K, reS p ec tively. The CC coherent pion cross section is assigned an error 
of 100% due to the fact that the CC coherent pion cross section had only 90% confidence 
upper limits for sub-GeV neutrino energies at the time of this analysis. In addition, when 
included in the MiniBooNE pion production fits, the data are consistent with the nominal 
NELIT model at lcr and with zero cross section at 2a. The NC coherent pion production 
data [76] differ from NEUT by 15%, within the measurement uncertainty of 20%. To account 
for the difference and the uncertainty, we conservatively assign a 30% overall uncertainty to 

rrNCcohir 

The anti-neutrino/neutrino cross section ratios are assigned an uncertainty of 40%. This 
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is a conservative estimate derived from doubling the maximum deviation between the energy- 
dependent MiniBooNE CCQE neutrino cross section and the RFG model assuming an axial 
mass of M% e = 1.03 GeV/c 2 , which was 20%. 

For v e CCQE interactions, there may be some effects that are not accounted for in the 
NEUT model, such as the existence of second class currents, as motivated in Ref. [77] . The 
dominant source of uncertainty is the vector component, which may be as large as 3% at 
the T2K beam peak, and thus is assigned as an additional error on v e CCQE interactions 
relative to CCQE interactions. 


Table [VIT] shows energy-dependent normalization parameters for CCQE and CC17T inter¬ 
actions which are included to account for possible discrepancies in the model as suggested, 
for example, by the difference between the MiniBooNE and NOMAD iza results. As men¬ 
tioned above, the uncertainties for x® E and xf c ' l7r are assigned from our study of MiniBooNE 
data. The remaining CCQE energy regions are assigned a 30% uncertainty to account for the 
aforementioned discrepancy while x ECl7r has a 40% uncertainty assigned since it is necessary 
to extrapolate from the MiniBooNE CCl7r + inclusive measurement at 2 GeV. 

The NCOther category consists of NCEL, NC resonant production where the resonance 
decays to a nucleon and kaon, eta, or photon, and NCDIS interactions. For fits to the ND280 
data and v e analyses at SK, resonant production that produces a nucleon and charged 
pion is also included in the NCOther definition, though kept separate in other analyses. 
NCOther interactions have a 30% normalization error assigned to them, which is given to 
the parameters x NCOth and x NCl7r± . 


D. Alternative Models 

As mentioned above, NEUT’s default model for CCQE assumes an RFG for the nuclear 
potential and momentum distribution of the nucleons. An alternative model, referred to as 
the “spectral function” (SF) [75], appears to be a better model when compared to electron 
scattering data. SF is a generic term for a function that describes the momentum and energy 
distributions of nucleons in a nucleus. In the model employed in [79], the SF consists of a 
mean-held term for single particles and a term for correlated pairs of nucleons, which leads to 
a long tail in the momentum and binding energy. It also includes the nuclear shell structure 
of oxygen, the main target nucleus in the T2K far detector. The difference between the 
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RFG and SF models is treated with an additional systematic parameter. 

At the time of this analysis, the SF model had not been implemented in NEUT, so the 
NuWro generator [SO] was used for generating SF interactions with the assumption that a 
NEUT implementation of SF would produce similar results. The SF and RFG distributions 
were produced by NuWro and NEUT, respectively, for and u e interactions on both carbon 
and oxygen, while using the same vector and axial form factors. 

The ratio of the SF and RFG cross sections in NuWro is the weight applied to each NEUT 
CCQE event, according to the true lepton momentum, angle, and neutrino energy of the 
interaction. Overall, this weighting would change the predicted total cross section by 10%. 
Since we already include in the oscillation analysis an uncertainty on the total CCQE cross 
section, the NuWro cross section is scaled so that at E u = 1 GeV it agrees with the NEUT 
CCQE cross section. 

A parameter xsf is included to allow the cross section model to be linearly adjusted 
between the extremes of the RFG (xsf= 0) and SF (x 5 j?=l) models. The nominal value 
for xsf is taken to be zero, and the prior distribution for xsf is assumed to be a standard 
Gaussian (mean zero and standard deviation one) but truncated outside the range [0,1]. 


E. Summary of cross section systematic parameters 


All the cross section parameters, x, are summarized in Tab. |VII[ including the errors 
prior to the analysis of near detector data. They are categorized as follows: 


1. Common between ND280 and SK; constrained by ND280 data. The parameters which 
are common with SK and well measured by ND280 are A4® E , M^ ES and some nor¬ 
malization parameters. 


2. Independent between ND280 and SK, therefore unconstrained by ND280 data. The 
parameters pf , Eb and SF are target nuclei dependent and so are independent between 
ND280 ( 12 C) and SK ( 16 0). 


3. Common between ND280 and SK, but for which ND280 data have negligible sensitivity, 

are 


so no constraint is taken from ND280 data. The remaining parameters in Tab. VII 


not expected to be measured well by ND280 and therefore are treated like independent 
parameters. 
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We define x n to be the set of cross section systematic parameters which are constrained by 
ND280 data (category 1), to distinguish them from the remaining parameters x s (categories 
2 and 3). 
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TABLE VII: Cross section parameters x for the ND280 constraint and for the SK 
oscillation fits, showing the applicable range of neutrino energy, nominal value, and prior 
error. The category of each parameter describes the relation between ND280 and SK and 
is defined in Sec. IIII El Parameters marked with an asterisk are not included in the 
parametrization for the appearance analysis. 


Parameter 

E v GeV Range 

units 

Nominal 

Error 

Category 

m q e 

all 

GeV/c 2 

1.21 

0.45 

1 

xf E 

0 <E V < 1.5 


1.0 

0.11 

1 

r QE 

1.5 <E V < 3.5 


1.0 

0.30 

1 

r QE 

X 3 

E v > 3.5 


1.0 

0.30 

1 

PF 12 C 

all 

MeV/c 

217 

30 


2 

E b 12 C * 

all 

MeV/c 

25 

9 


2 

PF 16 0 

all 

MeV/c 

225 

30 


2 

E b 16 0 * 

all 

MeV 

27 

9 


2 

xsf for C 

all 


0 (off) 

1 (on) 

2 

xsf for 0 

all 


0 (off) 

1 (on) 

2 

m r es 

all 

GeV/c 2 

1.41 

0.22 

1 

rpCC Ie 

X 1 

0 <E V < 2.5 


1.15 

0.32 

1 

rpCCln 

x 2 

E v > 2.5 


1.0 

0.40 

1 

x NC1tt° 

all 


0.96 

0.33 

1 

gCCcohn 

all 


1.0 

1.0 

3 

x CCOth 

all 


0.0 

0.40 

3 

j.N'ClV 

all 


1.0 

0.30 

3 

gNCcohn 

all 


1.0 

0.30 

3 

x NCOth 

all 


1.0 

0.30 

3 

W Shape 

all 

MeV/c 2 

87.7 

45. 

3 

3 

x PDD 

all 


1.0 

1.0 

3 

CC u e 

all 


1.0 

0.03 

3 

V jv 

all 


1.0 

0.40 

3 

x FSI 

all 


Section 

III B 1 


3 
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IV. NEAR DETECTORS 


Precision neutrino oscillation measurements require good understanding of the neutrino 
beam properties and of neutrino interactions. The two previous sections describe how we 
model these aspects for the T2K experiment and how we use external data to reduce model 
uncertainty. However, if only external data were used, the resulting systematic uncertainty 
would limit the precision for oscillation analyses. 

In order to reduce systematic uncertainty below the statistical uncertainty for the exper¬ 
iment, an underground hall was constructed 280 m downstream of the production target for 
near detectors to directly measure the neutrino beam properties and neutrino interactions. 
The hall contains the on-axis INGRID detector, a set of modules with sufficient target mass 
and transverse extent to continuously monitor the interaction rate, beam direction, and 
profile, and the off-axis ND280 detector, a sophisticated set of sub-detectors that measure 
neutrino interaction products in detail. 

This section describes the INGRID and ND280 detectors and the methods used to select 
high purity samples of neutrino interactions. The observed neutrino interaction rates and 
distributions are compared to the predictions using the beamline and interaction models, 
with nominal values for the systematic parameters. Section [V] describes how ND280 data 
are used to improve the systematic parameter estimates and compares the adjusted model 
predictions with the ND280 measurements. 

A. INGRID 

1. INGRID detector 

The main purpose of INGRID is to monitor the neutrino beam rate, profile, and center. 
In order to sufficiently cover the neutrino beam profile, INGRID is designed to sample 
the beam in a transverse section of lOmxlOm, with 14 identical modules arranged in two 
identical groups along the horizontal and vertical axes, as shown in Fig. [3j Each of the 
modules consists of nine iron target plates and eleven tracking scintillator planes, each made 
of two layers of scintillator bars (X and Y layers). They are surrounded by veto scintillator 
planes to reject charged particles coming from outside of the modules. Scintillation light 
from each bar is collected and transported to a photo-detector with a wavelength shifting 
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~10m 


FIG. 3: Overview of the INGRID viewed from beam upstream. Two separate modules are 
placed at off-axis positions off the main cross to monitor the asymmetry of the beam. 


fiber (WLS fiber) inserted in a hole through the center of the bar. The light is read out by 
a Multi-Pixel Photon Counter (MPPC) [8Tj attached to one end of the WLS fiber. A more 
detailed description can be found in Ref. [82] . 


2. Event selection 

Neutrino interactions within the INGRID modules are selected by first reconstructing 
tracks using the X and Y layers independently with an algorithm based on a cellular au¬ 
tomaton. Pairs of tracks in the X and Y layers with the same Z coordinates at the track 
ends are matched to form 3D tracks. The upstream edges of the 3D tracks in an event are 
compared to form a vertex. Events are rejected if the vertex is outside the fiducial volumes, 
the time is more than 100 ns from a beam pulse, or if there is a signal in the veto plane at 
the upstream position extrapolated from a track. 

This analysis [S3j significantly improves upon the original method established in 2010 [82]. 
The new track reconstruction algorithm has a higher track reconstruction efficiency and is 
less susceptible to MPPC dark noise. Event pileup, defined as more than one neutrino 
interaction occurring in a module in the same beam pulse, occurs in as many as 1.9% of 
events with interactions at the current beam intensity. The new algorithm handles pileup 
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events correctly as long as the vertices are distinguishable. For the full dataset, 4.81 x 10 6 
events are selected as candidate neutrino events in INGRID. The expected purity of the 
neutrino events in INGRID is 99.58%. 


3. Corrections 

Corrections for individual iron target masses and the background are applied in the same 
way as the previous INGRID analysis [82]. In addition, we apply corrections for dead 
channels and event pilcup which can cause events to be lost. There are 18 dead channels 
out of 8360 channels in the 14 standard modules and the correction factor for the dead 
channels is estimated from a Monte Carlo simulation. The correction factor for the event 
pileup is estimated as a linear function of the beam intensity, since the event-pileup effect is 
proportional to the beam intensity. The slope of the linear function is estimated from the 
beam data by combining events to simulate event pileup j83j. The inefficiency due to pilcup 
is less than 1% for all running periods. 


4- Systematic error 

Simulation and control samples are used to study potential sources of systematic error 
and to assign systematic uncertainties. The sources include target mass, MPPC dark noise 
and efficiency, event pilcup, beam-induced and cosmic background, and those associated 
with the event selection criteria. 

The total systematic error for the selection efficiency, calculated from the quadratic sum 
of all the systematic errors, is 0.91%. It corresponds to about a quarter of the 3.73% error 
from the previous analysis method [82]. The reduction of the systematic error results from 
the analysis being less sensitive to MPPC dark noise and event pileup, the improved track 
reconstruction efficiency, and more realistic evaluations of systematic errors which had been 
conservatively estimated in the previous analysis. 



FIG. 4: Daily event rate of the neutrino events normalized by protons on target. The error 
bars show the statistical errors. The horn current was reduced to 205 kA for part of Run 3. 


5. Results of the beam measurement 


Figure [4] shows the daily rates of the neutrino events normalized by POT. When the horn 
current was reduced to 205 kA due to a power supply problem, the on-axis neutrino flux 
decreased because the forward focusing of the charged pions by the horns becomes weaker. 
An increase by 2% and a decrease by 1% of event rate were observed between Runl and 
Run2, and during Run4, respectively. However, for all run periods with the horns operated 
at 250 kA, the neutrino event rate is found to be stable within 2% and the RMS/mean of 
the event rate is 0.7%. 


A Monte Carlo (MC) simulation that implements the beamline and neutrino interaction 
models described earlier, along with the INGRID detector simulation, is used to predict the 
neutrino event rate with the horns operating at 250 kA and 205 kA. The ratios of observed 
to predicted event rates, using the nominal values for the beamline and neutrino interaction 


systematic parameters, are: 

jydata 

^° c kA = 1.014 ± O.OOl(stat) ± 0.009(det syst), 

^250kA 

yydata 

Jg A = 1.026 ± 0.002(stat) ± 0.009(det syst), 

^205kA 


(3) 

(4) 


The uncertainties from the neutrino flux prediction and the neutrino interaction model are 
not included in the systematic errors. 

The profiles of the neutrino beam in the horizontal and vertical directions are measured 
using the number of neutrino events in the seven horizontal and seven vertical modules, 
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respectively. The observed horizontal and vertical profiles are fitted with separate Gaussian 
functions and the profile center is defined as the fitted peak positions. Finally, the neutrino 
beam direction is reconstructed as the direction from the proton beam target position to the 
measured profile center at INGRID using the result of accurate surveys of the proton beam 
target and the INGRID detectors. Figure [5] shows the history of the horizontal and vertical 
neutrino beam directions relative to the nominal directions as measured by INGRID and 
by the muon monitor. The measured neutrino beam directions are stable well within the 
physics requirement of 1 mrad. A 1 mrad change in angle changes the intensity and peak 
energy of an unoscillated neutrino beam at SK by 3% and 13 MeV, respectively. Because a 
misalignment in the proton beamline was adjusted in November 2010, the subsequent beam 
centers in the vertical direction are slightly shifted toward the center. A conservative esti¬ 
mate of the systematic error of the profile center is calculated by assuming that the detector 
systematic uncertainties for the neutrino event rate are not correlated between different 
INGRID modules. The average horizontal and vertical beam directions are measured as 

^beam _ q Q3Q o.Oll(stat) ± 0.095(det syst) mrad, (5) 

^beam _ q q;q -j_ o.012(stat) ± 0.105(det syst) mrad, (6) 

respectively. The neutrino flux uncertainty arising from possible incorrect modeling of the 
beam direction is evaluated from this result. This uncertainty, when evaluated without 
ND280 data, is significantly reduced compared to the previous analysis, as shown in Fig. [6j 
The horizontal and vertical beam width measurements are given by the standard de¬ 
viations of the Gaussians fit to the observed profiles. Figure [7] shows the history of the 

horizontal and vertical beam widths with the horns operating at 250 kA which are found 

to be stable within the statistical errors. The ratios of observed to predicted widths, using 
nominal values for the systematic parameters, are: 

1.015 ± O.OOl(stat) ± 0.010(det syst), (7) 

1.013 ± O.OOl(stat) ± 0.011(det syst), (8) 

for the horizontal and vertical direction, respectively. 
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FIG. 5: History of neutrino beam directions for horizontal (left) and vertical (right) 
directions as measured by INGRID and by the muon monitor (MUMON). The zero points 
of the vertical axes correspond to the nominal directions. The error bars show the 

statistical errors. 



FIG. 6: Fractional uncertainties of the flux at SK due to the beam direction uncertainty 
evaluated from the previous and this INGRID beam analyses. These evaluations do not 

include constraints from ND280. 


B. ND280 

In designing the experiment, it was recognized that detailed measurements of neutrino 
interactions near the production target and along the direction to the far detector would 
be necessary to reduce uncertainty in the models of the neutrino beam and of neutrino 
interactions. To achieve this, the T2K collaboration chose to use a combination of highly 
segmented scintillator targets and gaseous trackers in a magnetic spectrometer. Segmented 
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FIG. 7: History of neutrino beam width for horizontal (left) and vertical (right) directions 
for the horn 250 kA operation. The error bars show the statistical errors. 


active targets allow for the neutrino interaction to be localized and the trajectories of the 
charged particles to be reconstructed, and those passing through the gaseous trackers have 
their charge, momentum, and particle type measured. The targets and gaseous trackers 
are surrounded by a calorimeter to detect photons and assist in particle identification. The 
refurbished UAl/NOMAD magnet was acquired and its rectangular inner volume led to a 
design with rectangular sub-detectors. Spaces within the yoke allowed for the installation 
of outer muon detectors. 

The following sections describe the ND280 detector, its simulation, and the analyses used 
as input for the T2K oscillation analyses. 


1. ND280 detector 

The ND280 detector is illustrated in Fig. [8j where the coordinate convention is also 
indicated. The x and z axes are in the horizontal plane and the y axis is vertical. The origin 
is at the center of the magnet and the 0.2 T magnetic field is along the +£ direction. The 
z axis is the direction to the far detector projected onto the horizontal plane. 

The analyses presented in this paper use neutrino interactions within the ND280 tracker, 
composed of two fine-grained scintillator bar detectors (FGDs |B3j), used as the neutrino 
interaction target, sandwiched between three gaseous time projection chambers (TPCs [85]). 

The most upstream FGD (FGD1) primarily consists of polystyrene scintillator bars having 
a square cross section, 9.6 mm on a side, with layers oriented alternately in the x and y 
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FIG. 8: Sketch of the ND280 off-axis detector in an exploded view. A supporting basket 
holds the 7r° detector (POD) as well as the Time Projection Chambers (TPCs) and Fine 
Grained Detectors (FGDs) that make up the ND280 Tracker. Surrounding the basket is a 
calorimeter (ECal) and within the magnet yoke is the Side Muon Range Detector (SMRD). 


directions allowing projective tracking of charged particles. Most of the interactions in the 
first FGD are on carbon nuclei. The downstream FGD (FGD2) has a similar structure 
but the polystyrene bars are interleaved with water layers to allow for the measurement of 
neutrino interactions on water. The FGDs are thin enough that most of the penetrating 
particles produced in neutrino interactions, especially muons, pass through to the TPCs. 
Short-ranged particles such as recoil protons can be reconstructed in the FGDs, which 
have fine granularity so that individual particle tracks can be resolved and their directions 
measured. 

Each TPC consists of a field cage filled with Ar:CF 4 :iC 4 H 10 (95:3:2) inside a box filled with 
CO 2 . The +x and —x walls of the field cages are each instrumented with 12 MicroMEGAS 
modules arranged in two columns. The 336 mm x 353 mm active area for each MicroMEGAS 
is segmented into 1728 rectangular pads arranged in 48 rows and 36 columns, providing 3D 
reconstruction of charged particles that pass through the TPCs. The curvature due to the 
magnetic field provides measurements of particle momenta and charges and, when combined 
with ionization measurements, allows for particle identification (PID). 

The tracker is downstream of a n° detector (POD [86]) and all of these detectors are 
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surrounded by electromagnetic calorimeters (ECals [87]) and side muon range detectors 
(SMRDs [SB]). 

Data quality is assessed weekly. Over the entire running period, the ND280 data taking 
efficiency is 98.5%. For the analyses presented here, only data recorded with all detectors 
having good status are used, giving an overall efficiency of 91.5%. 


2. ND280 simulation 


A detailed simulation is used to interpret the data recorded by ND280. The neutrino 
flux model described in Sec. Ill Al is combined with the NEUT neutrino interaction model 


described in Sec. Ill A and a detailed material and geometrical description of the ND280 
detector including the magnet, to produce a simulated sample of neutrino interactions dis¬ 
tributed throughout the ND280 detector with the beam time structure. For studies of 
particles originating outside of the ND280 detector, separate samples are produced using a 
description of the concrete that forms the near detector hall and the surrounding sand. 

The passage of particles through materials and the ND280 detector response are modeled 
using the GEANT4 toolkit [H9]. To simulate the scintillator detectors, including the FGDs, 
we use custom models of the scintillator photon yield, photon propagation including reflec¬ 
tions and attenuation, and electronics response and noise [90]. The gaseous TPC detector 
simulation includes the gas ionization, transverse and longitudinal diffusion of the electrons, 
transport of the electrons to the readout plane through the magnetic and electric field, gas 
amplification, and a parametrization of the electronics response. 

Imperfections in the detector response simulation can cause the model to match the de¬ 
tector performance poorly, potentially generating a systematic bias in parameter estimates. 
After describing the methods to select neutrino interactions in the following section, we 
quantify the systematic uncertainty due to such effects with data/simulation comparisons 
in Sec. ITVB41 


3. ND280 zy x Tracker analysis 

We select an inclusive sample of z/ /x CC interactions in the ND280 detector in order to 
constrain parameters in our flux and cross section model. Our earlier oscillation analyses 
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divided the inclusive sample into two: CCQE-like and the remainder. New to this analysis 
is the division of the inclusive sample into three sub-samples, defined by the number of final 
state pions: zero (CC07r-like), one positive pion (CCl7r + -like), and any other combination 
of number and charge (CCOther-like). This division has enhanced ability to constrain the 
CCQE and resonant single pion cross section parameters, which, in turn, decreases the 
uncertainty they contribute to the oscillation analyses. 

The CC-inclusive selection uses the highest momentum negatively charged particle in an 
event as the yr candidate and it is required to start inside the FGD1 fiducial volume (FV) 
and enter the middle TPC (TPC2). The FV begins 58 mm inward from the boundaries 
of the FGD1 active volume in x and y and 21 mm inward from the upstream boundary 
of the FGD1 active volume in z, thereby excluding the first two upstream layers. The 
TPC requirement has the consequence of producing a sample with predominantly forward¬ 
going /i ~. Additional requirements are included to reduce background in which the start of 
the n~ candidate is incorrectly assigned inside the FGD1 FV, due to a failure to correctly 
reconstruct a particle passing through the FGD1 (through-going veto). The yr candidate is 
required to be consistent with a muon (muon PID requirement) based on a truncated mean 
of measurements of energy loss in the TPC gas [[85]. A similar PID has been developed 
for the FGD, which is not used for the muon selection, but is used in secondary particle 
identification [8T| . 

Events passing this selection comprise the CC-inclusive sample which is then divided 
into three exclusive sub-samples on the basis of secondary tracks from the event vertex. The 
names for these samples have the “-like” suffix to distinguish them from the corresponding 
topologies that are based on truth information. Those events with no additional TPC tracks 
consistent with being a pion or electron and with no additional FGD tracks consistent with 
being a pion, nor any time-delayed signal in the FGD which is consistent with a Michel 
electron, comprise the CCCD-like sample. Those events with one positive pion candidate 
in a TPC and no additional negative pions, electrons or positrons comprise the CCl7r + - 
like sample. The CCOther-like sample contains all other CC-inclusive events not in the 
CC07r-like or CCl7r + -like samples. 

In the simulation we find that the CC-inclusive sample is composed of 90.7% true CC 
interactions within the FGD fiducial volume, and 89.8% of the muon candidates are muons 


(the rest are mainly mis-identified negative pions). Table VIII shows the number of events 
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TABLE VIII: Number of events at each cut step, for data and for simulation (scaled to 

data POT) for the CC-inclusive sample. 


Requirement 

Data 

Simulation 

candidate starts within FGD1 FV and enters TPC2 

48731 

47752 

passes through-going veto 

34804 

36833 

passes muon PID requirement 

25917 

27082 


after each cut for data and simulation scaled to data POT, with systematic parameters set 
to their nominal values. 

Table |IX| shows that the CC07r-like sample is significantly enhanced in CCQE interactions, 
the CCl7r + -like sample in CC resonant pion interactions, and the CCOther-like sample in CC 
deep inelastic scattering (DIS) interactions. This division improves the constraints on several 
neutrino interaction model parameters. As shown in Tab. [X] the CCl7r + true topology is the 
most difficult to isolate. Most of the contamination in the CCl7r + -like sample comes from 
deep inelastic scattering events for which only one pion is detected and any other hadrons 
have escaped or have been lost to interactions in the surrounding material. 

Figures [9j [lOj [TlJ and 12 show the distributions of the muon momentum p /( and angle 
9n (with respect to the £-axis) for the CC-inclusive sample and each sub-sample. These are 
compared to the nominal simulation, broken down by true reaction type. 


4■ ND280 detector systematics 

In this section we explain how we use control samples to assess uncertainty in the modeling 
of FGD and TPC response and of neutrino interactions outside of the fiducial volume of the 
FGD. 

TPC systematic uncertainties are divided into three classes: selection efficiency, momen¬ 
tum resolution and PID. The efficiency systematic uncertainty arises in the modeling of the 
ionization, cluster finding (where a cluster is defined as a set of contiguous pads in a row or 
column with charge above threshold), track Ending, and charge assignment. This is assessed 
by looking for missed track components in control samples with particles that pass through 
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TABLE IX: Composition for the selected samples (CC-inclusive, CC07r-like, CCl7r + -likc, 
CCOther-like) according to the reaction types. 


True Reaction 

CC-inclusive 

CCOvr-like 

CClvr+-like 

CCOther-like 

CCQE 

44.6% 

63.3% 

5.3% 

3.9% 

Resonant pion production 

22.4% 

20.3% 

39.4% 

14.2% 

Deep inelastic scattering 

20.6% 

7.5% 

31.3% 

67.7% 

Coherent pion production 

2.9% 

1.4% 

10.6% 

1.4% 

NC 

3.1% 

1.9% 

4.7% 

6.8% 


0.5% 

0.2% 

1.7% 

0.9% 

v e 

0.3% 

0.2% 

0.4% 

0.9% 

Out of FGD1 FV 

5.4% 

5.2% 

6.6% 

4.1% 

Other 

0.05% 

0.03% 

0.04% 

0.2% 


TABLE X: Composition of the selected samples (CC-inclusive, CC07r-likc, CCl7r + -like, 
CCOther-like) divided into the true topology types. The non-z^ CC topology includes v e , 

z> M and NC interactions. 


True Topology 

CC-inclusive 

CC07r-like 

CCl7T+-like 

CCOther-like 

CCOtt 

51.5% 

72.4% 

6.4% 

5.8% 

CC17T + 

15.0% 

8.6% 

49.2% 

7.8% 

CCOther 

24.2% 

11.5% 

31.0% 

73.6% 

non-zy, CC 

4.1% 

2.3% 

6.8% 

8.7% 

Out of FGD1 FV 

5.2% 

5.2% 

6.6% 

4.1% 


all three TPCs. The single track-finding efficiency is determined to be (99.81 q 4%) for data 
and simulation for all angles, momenta and track lengths, and shows no dependence on the 
number of clusters for tracks with 16 clusters or more. The inefficiency due to the overlap 
from a second nearly collinear track is found to be negligible for both data and simulation, 
so this systematic uncertainty can be ignored. The same control samples are used to eval¬ 
uate the charge mis-identihcation systematic uncertainty. This systematic uncertainty is 
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Muon cos(0) 


FIG. 9: Muon momentum and angle distribution for the CC-inclusive sample. These are 
compared to the simulation, broken down into the different reaction types shown in 


Tab. IX and where non CC refers to NC, and u e interactions. All systematic 
parameters are set to their nominal values. 


evaluated by comparing data and simulation of the charge mis-identihcation probability as 
a function of momentum. This is found to be less than 1% for momenta less than 5 GeV/c. 

The momentum resolution is studied using particles crossing at least one FGD and two 
TPCs by evaluating the effect on the reconstructed momenta when the information from 
one of the TPCs is removed from the analysis. The inverse momentum resolution is found 
to be better in simulations than in data, typically by 30%, and this difference is not fully 
understood. A scaling of the difference between true and reconstructed inverse momentum 
is applied to the simulated data to account for this. Uncertainty in the overall magnetic 
held strength leads to an uncertainty on the momentum scale of 0.6%, which is confirmed 
using the range of cosmic ray particles that stop in the FGD. 

The TPC measurement of energy loss for PID is evaluated by studying high-purity control 
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FIG. 10: Muon momentum and angle distribution for the CC07r-like sample. These are 
compared to the simulation, broken down into the different reaction types, with all 
systematic parameters set to their nominal values. 


samples of electrons, muons and protons. The muon control sample has the highest statistics 
and is composed of particles from neutrino interactions outside the ND280 detector that pass 
through the entire tracker. For muons with momenta below 1 GeV/c, the agreement between 
data and simulation is good, while above 1 GeV/c the resolution is better in simulation than 
in data. Correction factors are applied to the simulation to take into account this effect. 

The performance for track finding in the FGD is studied separately for tracks which are 
connected to TPC tracks and tracks which are isolated in the FGD. The TPC-FGD matching 
efficiency is estimated from the fraction of through-going muons, in which the presence of 
a track in the TPC upstream and downstream of the FGD implies that a track should be 
seen there. The efficiency is found to be 99.9% for momentum above 200 MeV/c for both 
simulation and data. 

The FGD-only track efficiency is computed as a function of the direction of the track 
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FIG. 11: Muon momentum and angle distribution for the CCl7r + -like sample. These are 
compared to the simulation, broken down into the different reaction types, with all 
systematic parameters set to their nominal values. 


using a sample of stopping protons going from TPC1 to FGD1. This efficiency is found to 
be slightly better for data than simulation when cos 6^ < 0.9. A correction is applied to 
the simulation to account for this and the correction uncertainty is included in the overall 
detector uncertainty. 

The FGD PID performance is evaluated by comparing the energy deposited along the 
track with the expected energy deposit for a given particle type and reconstructed range 
in the FGD. We use control samples of muons and protons tagged by TPC1 and stopping 
in FGD1. The pull distributions (residual divided by standard error) for specific particle 
hypotheses (proton, muon or pion) for data and simulation are fitted with Gaussian distri¬ 
butions. To account for the differences in the means and widths of the distributions between 
data and simulation, corrections are applied to simulation and the correction uncertainty is 
included in the overall detector uncertainty. 
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FIG. 12: Muon momentum and angle distribution for the CCOther-like sample. These are 
compared to the simulation, broken down into the different reaction types, with all 
systematic parameters set to their nominal values. 


The Michel electron tagging efficiency is studied using a sample of cosmic rays that stop 
in FGD1 for which the delayed electron is detected. The Michel electron tagging efficiency is 
found to be (61.1 ± 1.9)% for simulation and (58.6 ± 0.4)% for data. A correction is applied 
to simulation and the correction uncertainty is included in the overall detector uncertainty. 

The uncertainty on the mass of the FGD, computed using the uncertainties in the size 
and density of the individual components, is 0.67% [84] . 

There is systematic uncertainty in the modeling of pion interactions traveling through the 
FGD. This is evaluated from differences between external pion interaction data jTT)l loTj and 
the underlying GEANT4 simulation. The external data do not cover the whole momentum 
range of T2K, so some extrapolation is necessary. Incorrect modeling can migrate events 
between the three sub-samples and for some ranges of momentum this produces the largest 
detector systematic uncertainty. 
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Ail out-of-fiducial volume (OOFV) systematic is calculated by studying nine different 
categories of events that contribute to this background. Examples of these categories are: 
a high energy neutron that creates a n~ inside the FGD that is mis-identified as a muon, 
a backwards-going n + from the barrel-ECal that is mis-reconstructed as a forward-going 
muon, and a through-going muon passing completely through the FGD and the TPC-FGD 
matching failed in such a way that mimics a FV event. Each of these categories is assigned a 
rate uncertainty (of 0 or 20%) and a reconstruction-related uncertainty. The reconstruction- 
related uncertainty is below 40% for all categories but one: we assign a reconstruction-related 
uncertainty of 150% to the high-angle tracks category, in which matching sometimes fails to 
include some hits that are outside the FGD FV. 

An analysis of the events originating from neutrino interactions outside the ND280 detec¬ 
tor (pit walls and surrounding sand) is performed using a dedicated simulation (sand muon 
simulation). The data/simulation discrepancy is about 10% and is included as a systematic 
uncertainty on the predicted number of sand muon events in the CC-inclusive sample. 

Pilcup corrections are applied to account for the inefficiency due to sand muons crossing 
the tracker volume in coincidence with a FV event. The correction is evaluated for each 
dataset separately and is always below 1.3%; the systematic uncertainty arising from this 
correction is always below 0.16%. 

Table |Xl| shows the full list of base detector systematic effects considered and the way each 
one is treated within the simulated samples to propagate the uncertainty. Normalization 
systematics are treated by a single weight applied to all events. Efficiency systematics are 
treated by applying a weight that depends on one or more observables. Finally, several 
systematics are treated by adjusting the observables and re-applying the selection. 

The base detector systematic effects are propagated using a vector of systematic param¬ 
eters d that scale the nominal expected numbers of events in bins of p ^-cos 6^ for the three 


selections, with the binning illustrated in Fig. 13 When a base systematic parameter is 
adjusted, di is the ratio of the modified to nominal expected number of events in bin i. The 
covariance of d due to the variation of each base systematic parameters is evaluated and 
the full covariance of d, Vd , is found by adding the individual covariances together. This 
covariance, and the observed number of events in the three samples in bins of p ^-cos 6^, 


shown in Fig. 13, are used by the subsequent analyses in order to constrain neutrino flux 
and interaction systematic parameters. 
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TABLE XI: List of base detector systematic effects and the way each one is treated within 
the simulated samples to propagate the uncertainty. Normalization systematics are treated 
with a signgle weight applied to all events. Efficiency systematics are treated by applying a 
weight that depends on one or more observables. Observable variation systematics are 
treated by adjusting the observables and re-applying the selection. 


Systematic effect 

treatment 

TPC tracking efficiency 

efficiency 

TPC charge misassignment 

efficiency 

TPC momentum resolution 

observable variation 

TPC momentum scale 

observable variation 

B Field distortion 

observable variation 

TPC PID 

observable variation 

TPC-FGD matching efficiency 

efficiency 

FGD tracking efficiency 

efficiency 

FGD PID 

observable variation 

Michel electron efficiency 

efficiency 

FGD mass 

normalization 

Pion secondary int. 

efficiency 

Out of Fiducial Volume 

efficiency 

Sand muon 

efficiency 

Pileup 

normalization 

TPC track quality requirements 

efficiency 
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FIG. 13: The p M -cos0 M binning for the systematic parameters d that propagate the base 
detector systematic effects are shown in the left figure for the three event selections. The 
binning for the observed number of events is shown in the right figure. For the CCl7r + -likc 
sample, the bin division at p M = 3.0 GeV/c is not used. 
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V. NEAR DETECTOR ANALYSIS 


In this section we explain how we use the large and detailed samples from ND280 in 
conjunction with models for the beam, neutrino interactions, and the ND280 detector to 
improve our predictions of the flux at SK and some cross section parameters. The systematic 
parameters for the beam model (b), binned in energy as shown in Fig. [lj the cross section 
model (x), listed in Tab. 


VII 


and detector model (d), illustrated in Fig. Il3l are used to 


describe the systematic uncertainties in the analysis. We use the three z/ /t CC samples 


described in Sec. IV B and external data discussed in Sec. Ill B and summarize our knowledge 
of the neutrino cross section parameters and unoscillated neutrino flux parameters with a 
covariance matrix, assuming that a multivariate Gaussian is an appropriate description. 


A. ND280 Likelihood 


The three ry, CC samples are binned in the kinematic variables p M and cos 9^, as shown 


in Fig. 13, and the observed and predicted number of events in the bins are used to define 
the likelihood, 


^bins 


C(h,x,S)=^p{Ni\Nf{b,x,S) 

i 

Nbins , x j\jd 


= c 


II 




(9) 


where Nf is the number of unoscillated MC predicted events and Nf is the number of data 
events in the ith bin of the CC samples, the second line assumes the Poisson distribution, 
and c is a constant. The number of MC predicted events, Nf(b,x,d), is a function of the 
underlying beam flux b, cross section x, and detector d parameters, and these parameters 
are constrained by external data as described in the previous sections. We model these 
constraints as multivariate Gaussian likelihood functions and use the product of the above 
defined likelihood and the constraining likelihood functions as the total likelihood for the near 
detector analysis. This total likelihood is maximized to estimate the systematic parameters 
and evaluate their covariance. In practice, the quantity —2 In C to t a i is minimized. Explicitly, 
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this quantity is: 


— 2 In Ctotai = constant + 


^bins 


( 10 ) 


2 V Nf(b,x,d) - N?\n{Nf(b,x,d)} 

2—1 

N b N b 

i= 1 i=i 
JV* N x 

i=l j=l 
N d N d 

+ 'L'L^- i ‘^"lbAd°i-d j ) 

i= 1 j=l 

where 6°, x°, and S' are the nominal values (best estimates prior to the ND280 analysis) 
and Vb, 14, and 14 are the covariance matrices of the beam, cross section, and detector 
systematic parameters. 


B. Fitting methods 

A reference Monte Carlo sample of ND280 events is generated using the models described 
in the previous sections and the nominal values for the systematic parameters. Predicted 
distributions for adjusted values of the systematic parameters are calculated by weighting 
each event of the Monte Carlo sample individually. For the flux parameters, the true energy 
and flavor of each MC event determine the normalization weight appropriate for that event. 
For the detector parameters, the reconstructed momentum and angle of the muon candidate 
are used. For cross section scaling parameters (e.g., xi), weights are applied according to 
the true interaction mode and true energy. For other cross section parameters (e.g., ), 

including the FSI parameters, the ratio of the adjusted cross section to the nominal cross 
section (calculated as a function of the true energy, interaction type, and lepton kinematics) 
is used to weight the event. The FSI parameters are constrained by a covariance matrix 


constructed by using representative points on the 1-cr surface for the parameters in Table IV 


The fit is performed by minimizing —2 In Ctotai using MINUIT [91]. Parameters not of 
interest to the oscillation analyses (e.g. ND280 detector systematic uncertainties) are treated 
as nuisance parameters. 
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FIG. 14: Comparison of the data and Monte Carlo distributions for muon momentum 
(top) and angle (bottom) in the CC07r-like sample, using the nominal and fitted values for 

the systematic parameters. 


Results 


The result of this analysis is a set of point estimates (g) and covariance (V g ) for the 
systematic scaling factors for the unoscillated neutrino flux at SK in bins of energy and 
flavor (b s ) and the cross section parameters which are constrained by ND280 data (x n ). 


Figures 14 [15} and 16 show the projected kinematic variable distributions of the three 
ND280 samples used in this analysis, comparing the data to the MC prediction for the two 
cases of using nominal values of the systematic parameters and using the best-fit values of 
the parameters. The MC distributions show better agreement with the data when using the 
best-fit values for the parameters, especially decreasing the prediction near the momentum 
peak and in the forward direction (cos 9^ close to 1). 

Figure [17] shows the values of the flux and cross section parameters that are constrained 


by the near detector analysis for the oscillation analyses; Table XII lists the flux parameters 


and Table XIII lists the values of the cross section parameters. These tables contain all of 
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Muon cos(0) 

FIG. 15: Comparison of the data and Monte Carlo distributions for muon momentum 
(top) and angle (bottom) in the CCl7r + -like sample, using the nominal and fitted values 

for the systematic parameters. 


the point estimates in g as well as the errors calculated as the square root of the diagonal 
of the covariance V g . One of the interesting features of the best-fit parameters is the dip in 
the flux parameters just below 1 GeV, which is near the peak of the T2K beam flux. This 
is particularly important, as this is the region of interest for oscillation analyses, and an 
incorrect prediction of the flux in this region can bias estimates of oscillation parameters. 
Another interesting point is the value of M EES , which is pulled to a much lower value than 
the external data constraint used in the fit. This highlights both the power of the ND280 
data, and the importance of the CCl7T + -like sample, which is dominant in determining this 
parameter. This selection is new to the ND280 analysis for the set of oscillation analyses 
reported in this paper, and provides an improved ability to use T2K data to constrain 
resonant interaction parameters. 

The predicted event rate at SK is given by the product of the flux, cross section, and 
detector efficiency, and the typical uncertainties of the flux and cross section parameters 
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FIG. 16: Comparison of the data and Monte Carlo distributions for muon momentum 
(top) and angle (bottom) in the CCOther-likc sample, using the nominal and fitted values 

for the systematic parameters. 


constrained by ND280 are 7-10%. The estimators of these flux and cross section parameters 
have a strong negative correlation, however, because they use the rate measurements in the 
near detector. As a result, their contribution to the SK event rate uncertainty is less than 
3%, significantly smaller than the individual flux and cross section parameter uncertainties. 

A cross-check to this analysis is performed by studying a selection of electron neutrino 
interactions in ND280 [92], and finds that the relative rate of selected electron neutrino 
events to that predicted by MC using the best-fit parameter values from this analysis is 
R(u e ) = 1.01 ±0.10. 
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FIG. 17: Prior and fitted values and uncertainties for the SK flux parameters (upper 
figure) and cross section parameters (lower figure) constrained by the near detector 
analysis for the oscillation analyses. Uncertainties are calculated as the square root of the 
diagonal of the relevant covariance matrix. The value of M® E and M EES are given in units 
of GeV/c 2 , and all other parameters are multiplicative corrections. 
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TABLE XII: Prior and fitted values and uncertainties for the near-detector-constrained SK 
flux parameters. All parameters are multiplicative corrections, and the uncertainties are 
calculated as the square root of the diagonal of the covariance matrix. 


Parameter 

Prior Value 

Fitted Value 

0.0-0.4 GeV 

1.004=0.12 

1.034=0.09 

Vn 0.4-0.5 GeV 

1.004=0.13 

1.024=0.09 

Vn 0.5-0.6 GeV 

1.004=0.12 

0.994=0.08 

z/ M 0.6-0.7 GeV 

1.00A0.11 

0.974=0.08 

z/„ 0.7-1.0 GeV 

1.004=0.13 

0.934=0.08 

Up 1.0-1.5 GeV 

1.004=0.12 

0.994=0.08 

1.5-2.5 GeV 

1.00=t0.10 

1.044=0.07 

z^ 2.5-3.5 GeV 

1.004=0.09 

1.054=0.06 

Vn 3.5-5.0 GeV 

l.OOAO.ll 

1.034=0.07 

z^ 5.0-7.0 GeV 

1.004=0.15 

0.984=0.07 

Vn >7.0 GeV 

1.004=0.19 

0.944=0.08 

z7 M 0.0-0.7 GeV 

1.004=0.13 

1.034=0.10 

0.7-1.0 GeV 

1.004=0.12 

1.014=0.09 

zq, 1.0-1.5 GeV 

1.004=0.12 

1.014=0.09 

1.5-2.5 GeV 

1.004=0.12 

1.034=0.10 

>2.5 GeV 

1.004=0.12 

1.014=0.11 

z/ e 0.0-0.5 GeV 

1.004=0.13 

1.034=0.10 

v e 0.5-0.7 GeV 

1.004=0.13 

1.014=0.09 

v e 0.7-0.8 GeV 

1.004=0.14 

0.984=0.11 

v e 0.8-1.5 GeV 

1.00A0.11 

1.004=0.07 

v e 1.5-2.5 GeV 

l.OOiO.lO 

1.024=0.07 

v e 2.5-4.0 GeV 

1.004=0.12 

1.004=0.07 

v e >4.0 GeV 

1.004=0.17 

0.954=0.08 

V e 0.0-2.5 GeV 

1.004=0.19 

1.014=0.18 

V e >2.5 GeV 

1.004=0.14 

0.964=0.08 
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TABLE XIII: Prior and fitted values and uncertainties for the near-detector-constrained 
cross section model parameters. The value of M® E and M EES are given in units of GeV/c 2 
and all other parameters are multiplicative corrections. The uncertainties are calculated as 
the square root of the diagonal of the covariance matrix. 


Parameter 

units 

Prior Value 

Fitted Value 

m% e 

GeV/c 2 

1.21T0.45 

1.24T0.07 

M% es 

GeV/c 2 

1.41T0.22 

0.96T0.07 

t Qe 

x i 


l.OOiO.ll 

0.97T0.08 

r QE 

x 2 


1.00T0.30 

0.93T0.10 

r QE 

x 3 


1.00T0.30 

0.85T0.11 

™CCl-7r 

x \ 


1.15T0.32 

1.26T0.16 

™CCl7T 

x 2 


1.00T0.40 

1.12T0.17 

x NCn ° 


0.96T0.33 

1.14T0.25 
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VI. FAR DETECTOR 


Precision measurements of neutrino oscillation by T2K rely on the capabilities of the far 
detector, most notably, its large target volume and acceptance and efficient discrimination 
between the primary leptons produced in z/ M and v e CC interactions. Additionally, since 
CCQE scattering interactions are expected to dominate at the energies below 1 GeV, accu¬ 
rate reconstruction of the parent neutrino energy is reliant upon accurate estimation of the 
lepton kinematics. Finally, the suppression of backgrounds, particularly those from NC and 
single-pion production processes, is needed. Here we discuss the performance of SK in this 
context, focusing on the event selections and the estimation of systematic uncertainties in 
the modeling of SK. 

Super-Kamiokande is a 50 kton water Cherenkov detector located in the Kamioka Obser¬ 
vatory, Gifu, Japan. It is divided into two concentric cylinders, an inner detector (ID) with 
11,129 inward-facing 20-inch photomultiplier tubes (PMTs) and an outer detector (OD), 
used primarily as a veto, which has 1885 outward-facing eight-inch PMTs. The ID PMTs 
view a 32 kton target volume and the OD collects light within a 2-m wide cylindrical shell 
surrounding the ID. The photocathode coverage of the ID is 40% and the space between 
PMTs is covered with a black plastic sheet to reduce reflection. To overcome its reduced 
photocathode coverage, reflective Tyvek® lines the inner and outer surfaces of the OD and 
each PMT is coupled to a 60 x 60 cm 2 wavelength-shifting plate to improve light collection. 

Cherenkov radiation from charged particles traversing the detector produces ring pat¬ 
terns recorded by the ID PMTs and is the primary tool for particle identification (PID). 
Due to their relatively large mass, muons passing through the detector are often unscat¬ 
tered and thereby produce clear ring patterns. Electrons, in contrast, scatter and produce 
electromagnetic showers, resulting in a diffuse ring edge. These differences in conjunc¬ 
tion with estimation of the Cherenkov opening angle enable efficient discrimination be¬ 
tween leptons. The probabilities to misidentify a single electron as a muon or a single 
muon as an electron are 0.7% and 0.8%, respectively, for typical lepton energies in T2K 
events. Since the recoil proton from CC interactions at T2K is usually below Cherenkov 
threshold, a single lepton is the dominant topology for beam-induced events at SK. For 
such isolated electrons (muons) the momentum and angular resolutions are estimated to be 
0.6% + 2.6%/^/P[GeV/c] (1.7% + 0.7%/-\/P[GeV/c]) and 3.0° (1.8°), respectively. Since the 
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start of T2K, SK has operated with upgraded electronics which provide lossless acquisition 
of all PMT hits above threshold. As a result the efficiency for tagging electrons from muon 
decays within the ID is 89.1%, an essential element of removing backgrounds containing 
sub-threshold muons or charged pions. Further details of the detector and its calibration 
may be found in [9j [93, f94] . 

Due to its large size, SK observes roughly ten atmospheric neutrino interactions per day 
within its fiducial volume. These neutrinos serve as control samples for the estimation of 
systematic errors. Similarly, although the detector is located at a depth of 2700 meters water 
equivalent, cosmic ray muons traverse the detector at approximately 3 Hz and together with 
their decay electrons provide an additional sample for systematic error evaluation. Details 
of these and other control samples are presented in the following subsections. 


A. Event Selection and Data Quality 

We define a sample of fully contained (FC) events whose Cherenkov light is deposited 
exclusively in the ID. PMTs in the OD that register light above threshold are referred to as 
“hit PMTs” and are grouped with neighboring hit PMTs to form clusters. If the largest such 
cluster contains more than 15 PMTs the event is rejected from the FC sample and included 
in the OD sample. Low energy (LE) events are removed by requiring that the total charge 
from the ID PMT hits in a 300 ns window be greater than 200 photoelectrons (p.e.), which 
corresponds to the charge observed from a 20 MeV electromagnetic shower. Events are also 
removed if a single ID PMT hit constitutes more than half of the total p.e. observed, in 
order to reject events due to noise. The final criterion rejects events that occur due to light 
from a discharge at the dynode of a PMT, known as “flasher” events. Such events have a 
broader timing distribution than neutrino interactions and tend to form a repeated pattern 
of light. A total of 18 events were rejected as flashers from all run periods, although from 
event timing information and visual scans we are confident that all are in fact due to beam 
neutrino interactions. Nevertheless, these events are discarded and the resulting selection 
inefficiency is taken into account. 

Events are timed with respect to the leading edge of the beam spill, taking into account 
the time of flight of the neutrino and myriad other sources of delay [9HI95]. Figure [18] shows 
the event timing (AT 0 ) distribution for all ID, OD, and LE events within ±500 /is of the 
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FIG. 18: ATo distribution of all FC, OD, and LE events within ±500 /is of the expected 
beam arrival time. The histograms are stacked in that order. 


beam arrival time. There is a clear peak near ATo = 0 for the FC sample. Eleven FC 
events have been observed outside the spill window. Using data collected with no beam we 
estimate the expected number of these events to be 5.85, mainly low energy events. ATo 
is corrected to take into account the neutrino interaction vertex position and the photon 
time-of-flight from the vertex to the PMTs. FC events within the spill window can be seen 


in Fig. 19 where the beam structure with eight bunches is clearly visible. The dotted lines 
represent the fitted bunch center times with a fixed bunch interval of 581 ns. For an event 
to be incorporated into the analysis, ATo must lie between —2 to 10 /is. 

A fiducial volume is defined within the ID, 2 m away from the detector wall, with a 
fiducial mass of 22.5 kton. Events whose vertex is reconstructed within this volume and 
with visible energy (T vis ) greater than 30 MeV are selected into the fully contained fiducial 
volume sample (FCFV). Visible energy is defined as the energy of an electromagnetic shower 
that produces the observed amount of Cherenkov light. We observe 377 events classified 
as FCFV. The expected number of background events from non-beam related sources in 
accidental coincidence is estimated to be 0.0085. 

Charged current interactions (v ± N —y l~ ± X ) in the narrow energy range of the T2K 
beam tend to produce single ring events at SK because most of the particles produced, 
except for the primary lepton, do not escape the nucleus or are below detection threshold. 
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FIG. 19: ATq distribution of all FC events within the beam spill window. 


The energy of the incoming neutrino can be calculated assuming the kinematics of a CCQE 
interaction and neglecting Fermi motion: 


rec m l - ( m *» - E b ) 2 - m 2 + 2 (m n - E b )Ei 
E,, = — 


( 11 ) 


2(m n - E b - Ei + pi cos 9 t ) 

where E™ c is the reconstructed neutrino energy, rn p is the proton mass, m n the neutron 
mass, mi the lepton mass and E b = 27MeV is the binding energy of a nucleon inside 16 0 
nuclei. E t , pi and 9i are the reconstructed lepton energy, momentum, and angle with respect 
to the beam, respectively. The selection criteria for both v e CC and v p CC events were fixed 
using MC studies before being applied to data. Events are determined to be e-like or //-like 
based on the PID of the brightest Cherenkov ring. The PID of each ring is determined by 
a likelihood incorporating information on the charge distribution and the opening angle of 
the Cherenkov cone. 


We select v e CC candidate events using the criteria listed in Tab. |XIV[ The E vis require¬ 
ment removes low energy NC interactions and electrons from the decay of unseen parents 
that are below Cherenkov threshold or fall outside the beam time window. The 7T°-like event 


rejection uses an independent reconstruction algorithm which is described in Sec. VIB We 
require E™ c < 1.25 GeV since above this energy the intrinsic beam v e background is dom¬ 
inant. The numbers of events remaining after successive selection criteria for a simulation 


sample produced with a nominal set of oscillation parameter values are shown in Tab. XIV 
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TABLE XIV: Event reduction for the v e CC selection at the far detector. The numbers of 
expected MC events divided into four categories are shown after each selection criterion is 
applied. The MC expectation is based upon three-neutrino oscillations for sin 2 20 2 3 = 1.0, 
= 2.4 x 10 _3 eV 2 /c 4 , sin 2 20 13 = 0.1, 5cp = 0 and normal mass hierarchy (parameters 
chosen without reference to the T2K data). 

(1) There is only one reconstructed Cherenkov ring 

(2) The ring is e-like 

(3) The visible energy, E v i s , is greater than 100 MeV 

(4) There is no reconstructed Michel electron 

(5) The reconstructed energy, E ^ ec , is less than 1.25 GeV 

(6) The event is not consistent with a 7r° hypothesis 



MC total 

^ ^ 

CC 

V e + Ve 

CC 

V + V 

NC 

u /t ->• v e 

CC 

interactions in FV 

656.83 

325.67 

15.97 

288.11 

27.07 

FCFV 

372.35 

247.75 

15.36 

83.02 

26.22 

(1) single ring 

198.44 

142.44 

9.82 

23.46 

22.72 

(2) electron-like 

54.17 

5.63 

9.74 

16.35 

22.45 

(3) E vis > 100 MeV 

49.36 

3.66 

9.68 

13.99 

22.04 

(4) no Michel election 

40.03 

0.69 

7.87 

11.84 

19.63 

(5) E r u ec < 1250 MeV 

31.76 

0.21 

3.73 

8.99 

18.82 

(6) not 7r°-like 

21.59 

0.07 

3.24 

0.96 

17.32 


After all cuts 28 events remain in the v e CC candidate sample. A Kolmogorov-Smirnov (KS) 
test of the accumulated events with accumulated POT is compatible with a constant rate 
with a p-value of 0.7. 


We select CC candidate events using the selection criteria shown in Tab. XV The 
momentum cut rejects charged pions and misidentihed electrons from the decay of unob¬ 
served muons and pions. We require fewer than two Michel electrons to reject events with 
additional unseen muons or pions. After all cuts are applied, 120 events remain in the v^ 
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TABLE XV: Event reduction for the CC selection at the far detector. The numbers of 
expected MC events divided into four categories are shown after each selection criterion is 
applied. The MC expectation is based upon three-neutrino oscillations for sin 2 26*23 = 1-0, 
Am 2 2 = 2.4 x 10~ 3 eV 2 /c 4 and normal mass hierarchy (parameters chosen without reference 

to the T2K data). 

(1) There is only one reconstructed Cherenkov ring 

(2) The ring is /r-like 

(3) The reconstructed momentum, p^, is greater than 200 MeV/c 

(4) There are less than two reconstructed Michel electrons 



MC total 

CCQE 

+ I'm 

CC nonQE 

u e + Z7 e 

CC 

V + V 

NC 

interactions in FV 

656.83 

111.71 

213.96 

43.05 

288.11 

FCFV 

372.35 

85.55 

162.20 

41.58 

83.02 

(1) single ring 

198.44 

80.57 

61.87 

32.54 

23.46 

(2) muon-like 

144.28 

79.01 

57.80 

0.35 

7.11 

(3) > 200 MeV/c 

143.99 

78.84 

57.77 

0.35 

7.04 

(4) Xjyiichei—e A 1 

125.85 

77.93 

40.78 

0.35 

6.78 


CC candidate sample. 


Figure 20 shows the candidate event spectra for the appearance {y e ) and disappearance 
(zyj channels. We monitor the vertex distributions of the candidate event samples for signs of 
bias that might suggest background contamination. Figure [21] shows the vertex distribution 
of the v e CC candidate events in the SK tank coordinate system. We observe no unexpected 
clustering and combined KS tests for uniformity in r 2 and z yields a p-value of 0.6. 


B. 7T° Rejection with the New Event Reconstruction Algorithm 


As mentioned in the previous section, in order to select v e CC events, we require that 


only one electron-like ring is reconstructed. The u e CC selection criteria 1-5 in Tab. XIV 
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(a) (b) 

FIG. 20: The reconstructed energy spectra of the observed v e (a) and (b) CC candidate 
events assuming CCQE interaction kinematics. The data are shown as points with 
statistical error bars and the shaded, stacked histograms are the MC predictions, and the 
rightmost bin includes overflow. The expectation is based on the following oscillation 
parameters: sin 2 26*13 = 0.1, sin 2 20 2 3 = 1.0, Scp = 0, Am | 3 = 2.4 x 10 - 3 eV 2 /c 4 and normal 

mass hierarchy. 


are based on the information provided by SK event reconstruction software which has been 
used at SK for atmospheric neutrino and nucleon decay analyses[1] and, as shown in the 
table, we reject most of the background events by these selection cuts. The v e appearance 
signal purity is 59.3% and the selection efficiency for the signal is 71.8%. The remaining 
backgrounds are predominantly NC single 7 r° events, as one of the two decay ys from a 7 r° 
is occasionally missed and the other 7 forms an electron-like ring. 


In order to reject such 7 r° events, we employ a new event reconstruction algorithm which is 
based on the methods developed by MiniBooNE [96j. The new algorithm adopts a maximum 
likelihood method to reconstruct particle kinematics in the SK detector. For a given event, 
we construct a likelihood function which uses the observed charge and time information from 
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(a) 



Vertex R 2 (cm 2 ) 

(b) 


FIG. 21: Two-dimensional vertex distributions of the observed v e CC candidate events in 
{x, y ) and (r 2 = x 2 + y 2 , z). The arrow indicates the neutrino beam direction and the 
dashed (blue online) line indicates the fiducial volume boundary. Events indicated by open 
square markers passed all of the v e selection cuts except for the fiducial volume cut. 


the PMTs: 

unhit 

l[x )=n Pj (unhit | x) 

j ( 12 ) 

hit v ' 

- Pi (unhit \x)}f q (qi\x)ft(ti\x). 

i 

In the equation, x represents particle track parameters such as the vertex, direction, and 
momentum which are to be estimated. The Erst index j runs over the PMTs which do not 
register a hit, and for each of such PMTs the conditional probability Pj (unhit|cc) of not 
registering a hit given x is evaluated. For each PMT which does register a hit, in addition 
to the hit probability, we calculate the probability density f q (qi\x) of observing charge qi 
as well as the probability density f t (ti\x) of the hit occurring at time fj. The estimated 
track parameters, x , are those that maximize the likelihood function. For every event 
we construct and maximize the likelihood assuming several different particle hypotheses, 
and particle identification is done using ratios of the maximum likelihoods for the different 
hypotheses. 

In this analysis, we use a single electron hypothesis and a i r° hypothesis for 7T° rejection. 
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FIG. 22: 2D distributions of the logarithm of the likelihood ratio ln(L„.o/L e ) vs. the 
reconstructed invariant mass m 77 , for signal v e CCQE(left) and background NC 7 T°(right) 
events. The diagonal line indicates the 7 r° rejection criterion, and events lying above the 
line are rejected as 7T° background. The size of each box is proportional to the number of 
events the bin. The two figures use the same scale for representing the number of events 

and are normalized to the same POT. 


The single electron hypothesis has seven parameters which are the initial vertex position, 
time, direction, and momentum. Since a 7 r° decays into two 7 s and produces two electron¬ 
like Cherenkov rings, the 7 T° hypothesis is constructed by combining the charge and time 
contributions from two electron tracks which point back to a common vertex. In addition 
to the common vertex and the directions and momenta of the two 7 tracks, each track has 
an additional free parameter which shifts its origin along its direction in order to account 
for photon conversion points. The 7 T° hypothesis therefore has twelve parameters. 

In order to distinguish signal v e CC events from 7 T° background events, we use the maxi¬ 
mum likelihood values of the electron hypothesis L e and the 7 r° hypothesis L n 0 as well as the 
reconstructed invariant mass m 77 obtained from the 7 r° hypothesis. Figure [22] shows the 2D 
distributions of the logarithm of the likelihood ratio ln(L 7 r o/L e ) vs. m 77 for signal v e CCQE 
and background NC 7 r° events which satisfy the v e selection criteria 1-5, produced by MC. 
We see a clear separation between the two event types, and we accept an event as a v e CC 
candidate if it satisfies ln(L 7 r o/L e ) < 175 — 0.875 x m 77 [MeV/c 2 ], which is indicated by the 


diagonal line in the plots. As shown in Tab. XIV the remaining NC background is reduced 
by roughly a factor of nine by introducing the 7 r° rejection cut. After the cut, the purity 
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FIG. 23: Efficiencies for rejecting NC 7 r° events for the previous and the new 7 r° rejection 
methods, plotted in bins of the energy of the less energetic 7 . 


and the selection efficiency for the v e appearance signal are 80.2% and 66.1% respectively. 

In earlier published T2K v e appearance analysis results PEZ], we used a 7 r° rejection 
method which is different from what is described above [98]. To demonstrate the improve¬ 


ment over the previous method, Fig. 23 shows the efficiency for rejecting NC 7T° events for 
the two methods, plotted as a function of the energy of the less energetic 7 . In calculating 
the efficiencies, only the events which satisfy the v e selection criteria 1-5 are included. As the 
figure indicates, the rejection efficiency by the new method remains high even in cases where 
the energy of one of the two 7 s is low. By employing the new method, we have reduced the 
7 T° background remaining in the final u e CC candidate event sample by 69% relative to the 
previous method. 


C. Systematic uncertainty 

This section describes the studies and treatment of uncertainty in modeling the SK detec¬ 
tor that lead to systematic uncertainty in estimating the selection efficiency and background 
for the oscillation samples. We use SKDETSIM PE], a GEANT3-derived simulation of the 
SK detector, to model the propagation of particles produced by neutrino interactions. The 
GCALOR physics package is used to simulate hadronic interactions in water owing to its 
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ability to reproduce pion interaction data around lGeV/c. However for pions with momen¬ 
tum below 500MeV/c, custom routines are employed based on the cascade model used by 
NEUT to simulate interactions of final state hadrons. SKDETSIM incorporates the propa¬ 
gation of photons in water, subject to absorption, Rayleigh scattering, and Mie scattering. 
The simulation of these processes is tuned using laser calibration sources in situ [93] . 

Control samples that are not related to the T2K beam spills are used to assess systematic 
uncertainty, including muons and neutrinos produced from cosmic ray interactions with the 
atmosphere (cosmic ray muons and atmospheric neutrinos) and combinations of simulated 
and cosmic ray data (hybrid-7T° sample). As described below, cosmic ray muons are used to 
evaluate the systematic uncertainty due to the fully-contained (FC), fiducial-volume, and 
decay-electron requirements. Atmospheric neutrinos are used to assess uncertainty from the 
ring counting, particle identification, and n° rejection. The hybrid-7T° sample is used to 
study the SK response to 7r 0, s. The uncertainties due to energy scale, modeling of pion final 
state interactions (FSI) and secondary interactions (SI) are evaluated separately. 

Cosmic ray muon samples are used to estimate uncertainties related to the FC, fiducial- 
volume and decay-electron requirements, for the selections of both v e and CC candidates. 
The error from the initial FC event selection is 1% and is dominated by the event-by-event 
flasher rejection cut. The uncertainty in the fiducial volume is estimated to be 1% using 
the vertex distribution of cosmic ray muons which have been independently determined to 
have stopped inside the ID. The uncertainty due to the Michel electron tagging efficiency is 
estimated by comparing cosmic ray stopped-muon data and MC. This uncertainty is applied 
based on the fraction of events with true Michel electrons in the T2K beam MC. The rate of 
falsely identified Michel electrons is estimated from MC and 100% uncertainty in that rate 
is assumed. Overall, the event rate uncertainty related to the decay-electron requirements 
is small. For the v e CC candidate sample, it is 0.2% for v e CC events and 0.4% for CC 
and NC events. For the CC candidate sample it is 1.0%. 

Other studies of systematic uncertainty in SK modeling divide simulated events into 
categories according to their final state (FS) topologies, with the criteria shown in Tab. 

These topologies do not correspond exactly with true interaction modes due to subsequent 
interactions within the nucleus or with neighboring nuclei or because one or more particles 
are produced below Cherenkov threshold. 

Atmospheric neutrino data are used to assess possible mismodeling of the ring count- 


XVI 
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TABLE XVI: Criteria for categorization of simulated events by final state topology for 
systematic studies. N x is the number of particles of type x and the number of charged 
pions (N n ±) and protons (Np) only includes those particles produced with momentum 
above Cherenkov threshold set at 156.0 MeV/c and 1051.0 MeV/c respectively. 


Event type 
CC le 
CC e other 
CC 1/j 


MC truth selection criteria 

v e CC and N n o = 0 and N n ± = 0 and Np = 0 

v e CC and not v e CCle 

v fL CC and N n o = 0 and N w ± = 0 and Np = 0 


CC n other 
CC n 7T° other 
NC It r° 

NC 7 T° other 


Vfx CC and N n o = 0 
jy CC and N n o > 0 

NC and not NC ly and N n o = 1 and N n ± = 0 and Np = 0 
NC and not NC I 7 and N n 0 > 1 and not NC l 7 r° 


NC I 7 NEUT truth 

NC I 7 r 1 * 3 NC and not NC I 7 and N n 0 = 0 and N n ± = 1 and Np = 0 

NC other NC and not NC I 7 and not NC l 7 r° and not NC Itt^ and not NC 7 r° other 


ing (RC), particle identification, and 7 r° rejection for the first four FS topologies shown in 


Tab. XVI Atmospheric neutrino samples fully contained within the fiducial volume and with 
E v i S > 30 MeV are divided into CCQE and CC non-quasi-elastic (CC 11 QE) enriched samples 
using the number of Michel electrons and the visible energy. These samples are further split 
into “core” samples of events which pass all of the requirements and tail samples of events 
which fail only one requirement. An additional background sample is included, enhanced in 


NC 7 T°. These samples, 13 in total, are summarized in Tab. XVII and are binned in E vis , for 
E vis < 30 GeV. 

In order to adjust the modeling of ring counting, particle identification, and 7 T° rejection, 
a set of parameters is defined to alter the cut values for these three classifiers. Separate pa¬ 


rameters are used for the first four FS topologies in Tab. XVI and for each visible energy bin 
within those topologies. By adjusting these parameters, simulated events, generated accord¬ 


ing to models of the atmospheric neutrino flux, migrate between the branches in Tab. XVII 


thus changing the efficiency for true CC le and CC e other (CC 1/j and CC /j other) events 
in the v e (z/^) core samples. Using the observed numbers of core and tail data events in each 
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TABLE XVII: SK atmospheric neutrino control samples. The parent sample is defined to 
be fully contained and in the fiducial volume. This parent sample is divided into four sets 
of core and tail samples and one background (BG) control sample. The main difference 
between CCQE and CCnQE is the number of decay-e N dcy _ e cut, which is based on the hit 
time distribution. The distance from the expected muon stopping point to the nearest 
decay-e, D dcy _ e , is used to select high purity n /( CCQE and CCnQE samples. The BG 
sample is enriched in NC 7r° to constrain the NC normalization. 


Type of control sample 


Branch of control sample 

sample 

RC cut 

PID cut 

/ D dcy _ e cut 



Core 

1R 

& e-like 

& not 7r°-like 

u e CCQE 

Ndcy—e — 0 

RC tail 

> 1R 

& e-like 

& not 7r°-like 

enriched 

& E V i S > 100 

PID tail 

1R 

& fi- like 

& not 7r°-like 



7r° tail 

1R 

& e-like 

& 7r°-like 



Core 

1R 

& e-like 

& not 7r°-like 

u e CCnQE 

Ndcy—e — 1 

RC tail 

> 1R 

& e-like 

& not 7r°-like 

enriched 

& E V i S > 100 

PID tail 

1R 

& p- like 

& not 7r°-like 



7T° tail 

1R 

& e-like 

& 7r°-like 

uy CCQE 

Ndcy—e 1 

Core 

1R 

& /I- like 

& D dcy _ e < 80 cm 
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II 
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visible energy bin, a likelihood function is defined and marginalized over the neutrino flux, 
neutrino interaction systematic parameters and cut adjustment parameters, using a Markov 
Chain Monte Carlo. The marginalized likelihood is used to estimate corrected efficiencies 
for the four FS topologies in bins of E v - ls and their covariance. The observed differences 
between the nominal and corrected efficiencies may indicate mismodcling of the detector 
response, so additional covariance is included, with the diagonal elements being the square 
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of these differences and the off-diagonal terms calculated by assuming full correlation. The 


correlations between the estimated efficiencies are shown in Fig. 24 

To evaluate the systematic uncertainty in modeling 7 r 0, s in SK, we construct a set of 
“hybrid- 7 T°” control samples. These events are constructed by overlaying one electron-like 
ring from the SK atmospheric neutrino or cosmic ray muon samples with one simulated 
photon ring. The simulated photon ring kinematics are chosen such that the momenta and 
opening angle between the two rings follow the decay kinematics of NC 7T° events from 
the T2K MC. Hybrid-7T° MC samples with both rings from the SK MC are produced to 
compare with the hybrid-7T° data samples and the difference in the fractions that pass the 
v e selection criteria is used to assign the systematic error. The difference could be due to 
incorrect modeling of scattered or reflected light from the higher energy ring which obscures 
the lower energy ring. In order to investigate this, we compare hybrid- 7 T° samples in which 
the electron constitutes the higher energy ring from the 7 r° decay with hybrid- 7 r° samples in 
which it constitutes the lower energy ring. For events with additional particles in the final 
state, we add a MC ring to the existing hybrid- n° samples and assume the dominant source 
of error comes from the detection of the lower energy photon from the 7T° decay. Relative 
uncertainties on the efficiency, calculated for 17 bins in reconstructed electron momentum 
and angle, are in the range 2-60%. Relative statistical errors, in the range 15-50%, are 
applied assuming no correlation between bins. 

Neutral-current interactions can produce a final state containing just a single photon via 
radiative decays of A resonances (NC ly). This is a background in the v e CC candidate 
sample because photons and electrons produce very similar charge patterns in the SK detec¬ 
tor. The uncertainty in the efficiency of selecting NC I 7 events is determined by comparing 
the efficiency of a single photon MC sample with that of a single electron MC. The difference 
is no more than 1 %. This error is added in quadrature to the uncertainty for the CC le FS 
topology described above to give the total uncertainty on the NC I 7 background in the v e 
CC candidate sample. 

Muon decay-in-flight events make up a small background in the v e CC candidate sample. 
Such events can be misidentihed because the decay electron is boosted in the direction of 
the parent muon and thus their Cherenkov rings can overlap. MC studies indicate that 
such events make up 19% of the background from interactions and its rate is assigned a 
16% selection uncertainty. The remaining background from v^ interactions are assigned a 
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FIG. 24: The correlations between the estimated efficiencies for the final state topologies 
CC le and CC e other for the v e CC event selection and CC 1/i and CC /./ other for the 
CC event selection. The upper figure shows the combinations with positive correlation, 
and the lower with negative correlation. The diagonal correlations (correlation = 1) are 

not shown. 
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conservative 150% error. A conservative 100% nncorrelated error is assigned to the NC 
and the NC other. 


For the CC candidate sample, the dominant NC backgrounds are NC and events 
with just a single proton (NC other). The relative uncertainty in this background, due to 
systematic uncertainties in ring counting and particle identification, is found to be 59%. The 
background from u e interactions in the CC candidate sample is assigned a conservative 
100 % error. 

All aspects of SK detector simulation that can affect the modeling of the SK candidate 
event selection described above are propagated using a vector of systematic parameters, 
s, which scale the nominal expected number of events in bins of the observable kinematic 
variables E rec or p — 6 for the true v interaction mode categories. The binning is shown in 


Tab. |XVIII| The covariance of these parameters, V s , is used to propagate the uncertainties 
in the detector simulation to the oscillation analyses. 

The energy scale uncertainty is estimated by comparing data with simulated samples 
spanning the momentum range 30MeV/c to 6GeV/c. Starting at the lowest energy, we use 
the reconstructed momentum spectrum of electrons produced by the decay of cosmic ray 
muons, the reconstructed mass of neutral pions from atmospheric neutrino interactions and 
cosmic ray muons that stop within the SK tank. The final uncertainty is 2.4%, independent 
of E v . 


Systematic uncertainties in pion interactions in the target nucleus (FSI uncertainties) 
and SK detector (SI uncertainties) are evaluated by varying pion interaction probabilities 
in the NEUT cascade model. In the NEUT sample we store the information necessary to 
recompute the pion cascade using modified interaction probabilities to weight each event. 
Altered CC sample distributions are produced using 16 representative points x FSI on the 
1 -er surface for the parameters. The covariance matrix V, which describes the variations in 
the number of events in the binned observables (A% due to the variation in x FSI , is given 
by 

1 16 

v„ = - £(JV ( (Sf") - M) - JVj) (13) 

U k=1 


The binning of this matrix is chosen to match that of the detector error covariance matrix 
shown in Table IXVIIII 

A simulation of photo-nuclear (PN) interactions is incorporated into the SK MC. The 





TABLE XVIII: Binning for the vector of SK detector systematic parameters s. Two 
schemes are defined since analyses use either E rec or p — 6 binning for the v e appearance 

channel. 


u e appearance: E rec ( GeV) 

Osc.r'eCC 

0-0.35-0.8-1.25 (3 bins) 

^cc 

0-0.35-0.8-1.25 (3 bins) 

v e CC 

0-0.35-0.8-1.25 (3 bins) 

NC 

0-0.35-0.8-1.25 (3 bins) 

u /t disappearance: E rec (GeV) 

z^CCQE 

0-0.4-1.1-30.0 (3 bins) 

n /t CCOther 0-30.0 (1 bin) 

v e 

0-30.0 (1 bin) 

NC 

0-30.0 (1 bin) 


v e appearance 

p (GeV/c) 9 (degree) 

| 

. 0-0.3 0-40-60-80-100-120-140-180 (7 bins) 

Osc.i/ e CC:i/^CC:i/ e CC:NC < 

I 0.3-0.7 0-40-60-80-180 (4 bins) 

0.7- 0-40-180 (2 bins) 

l/n disappearance 

E rec (GeV) 

i/^CCQE 

0-0.4-1.1-30.0 (3 bins) 

n /t CCOther 

0-30.0 (1 bin) 

V e 

0-30.0 (1 bin) 

NC 

0-30.0 (1 bin) 


model allows for the absorption of photons based on the measured cross section and assumes 
that there is no subsequent emission above Cherenkov threshold. A systematic uncertainty 
of 100% is assumed for the normalization of the PN cross section. 
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VII. OSCILLATION MODEL AND PARAMETER ESTIMATION 


The previous sections have described the T2K experiment and the way we model all 
elements of the experiment and neutrino interactions which are necessary to interpret our 
data, and how we use internal and external data to improve our models. In this section, 
we turn our attention to general aspects of estimating neutrino oscillation parameters from 
our data. The oscillation model is given along with the predictions for the probability for 
muon neutrino disappearance and electron neutrino appearance, the key observables for our 
experiment. We explain how we use external data for some of the oscillation parameters and 
the general approaches we use to estimate the remaining parameters. Finally, we characterize 


the importance of the different sources of systematic uncertainty. Sections VIII X describe 
the individual analyses and their results in detail. 


A. Oscillation model 

The Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix, U , defines the mixture of the 
mass eigenstates (v \, z/ 2 , and z/ 3 ) that make up each flavor state: 


(14) 


and it has become standard to parametrize this matrix, ignoring the Majorana phases, as: 
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(15) 


U = 


where Sij = sin 6 iv Cij = cos 9ij, and 5 = 5cp is the CP-violating phase. 

The ^-survival probability for a neutrino with energy E traveling a distance L is: 

-A z/ M ) = 

1-4 (s 2 12 c% 3 + sf 3 s% 3 cf 2 + 2si 2 si 3 s 23 ci 2 c 23 cos 5) s% 3 cj 3 sin 2 </> 31 
-4 (c 2 2 c 2 3 + sj 3 sl 3 sl 2 - 2si 2 Si 3 s 23 Ci 2 c 23 cos 5) s 2 23 cl 3 sin 2 0 32 

-4 (sj 2 c% 3 + s 2 13 s 2 23 c 2 2 + 2si 2 Si 3 s 23 Ci 2 c 23 cos (5) (c 2 2 c 2 3 + S? 3 s^ 3 sf 2 - 2 si 2 Si 3 s 23 Ci 2 c 23 cos (5) sin 2 (p 2 \ 

(16) 
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where 


'*7 


A mjjL 
4 E 


(17) 


in natural units and Amf - = rrif — m 2 is the difference in the squares of masses of eigenstates. 

The zvappearance probability, to first order approximation in matter effects, can be 
written as: 

P(Vl* -t Ve) = 

4c 2 13 s 2 13 s 2 23 sin 2 (j) 31 (l + a%(1 - 24 s)) 

+ 8 c? 3 Si 2 Si 3 s 2 3 (ci 2 c 23 cos 5 - s 12 s 13 s 23 ) cos 0 23 sin 0 3 1 sin 0 2i 

- 8 c 33 ci 2 c 23 si 2 si 3 s 23 sin <5 sin 0 32 sin 0 3 i sin 0 2i 

+ 4s 12 c 13 ( c 12 c 23 + S 12 S 23 S 13 “ 2ci 2 C 23 Si 2 S 23 Si 3 COS 5) sill 2 (j) 21 

-8c 2 13 s 2 13 sl 3 (1 - 24 ) cos 032 sin 031 

The effect on oscillation due to the density, p, of matter through which the neutrinos travel 
is included with the terms, a[eV 2 /c 4 ] = 7.56 x 10~ 5 p[g/cm 3 ]E^[GeV]. The corresponding 
I7 e -appearance probability is calculated by changing the sign of a and 5cp- Our analyses 
use the complete formulas, without approximating matter effects, to compute the oscillation 
probabilities. 

Since the neutrino mass hierarchy (MH) is not yet known, we parametrize the large mass 
splitting by |Am 2 | = A m 32 for normal hierarchy (NH, where m 3 is the largest mass) and 
|Am 2 1 = Am 2 3 for inverted hierarchy (IH, where m 3 is the smallest mass). 

It is not possible to estimate all of the oscillation parameters using only our measurements 
of z^-disappearance and z/ e -appearance. Instead, we estimate the four oscillation parameters, 
|Am 2 |, sin 2 0 23 , sin 2 0i 3 , 5cp, and the mass hierarchy, and use external measurements for the 
solar oscillation parameters, sin 2 di 2 and A m 21 , as we have negligible sensitivity to those. 
Figure [25] illustrates how our key observables depend on the two parameters, sin 2 0 23 and 5cp, 
for the two mass hierarchies. In this figure the neutrino energy is at the oscillation maximum 
(0.6 GeV), and the other oscillation parameters are fixed (solar parameters as established 
and sin 2 d 13 = 0.0243). To a good approximation, with our current dataset, 


in Sec. 


VIIB 


^-disappearance can be treated on its own to estimate 9 23 . The oscillation parameter 
dependence on zz e -appearance cannot be factorized, however. In order to estimate the full 
set of oscillation parameters and properly account for all uncertainties, it is necessary to do 
a joint analysis of z/^-disappearance and zz e -appearance. 
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FIG. 25: The P{v^ —* v^) survival probability and P{y M —> u e ) appearance probability for 
different values of sin 2 #23 and for Scp in the interval [— tt, n\ for normal (solid) and inverted 
(dashed) mass hierarchy. The highlighted dot on each ellipse is the point for Scp = 0 and 
Scp increases clockwise (anti-clockwise) for normal (inverted) mass hierarchy. The other 


oscillation parameter values are fixed (solar parameters as established in Sec. VIIB and 
sin 2 0i3 = 0.0243) and the neutrino energy is fixed to 0.6 GeV. 


B. External input for oscillation parameters 


Since our experiment is insensitive to the solar oscillation parameters, we fix them to the 
values sin 2 # 12 = 0.306 and = 7.5 x 10~ 5 eV 2 /c 4 from [§§] . As a check, the Bayesian 

analysis presented in Sec. [X] applies Gaussian priors with standard deviations (0.017 and 
0.2 x 10 _5 eV 2 /c 4 ) and finds that the uncertainties in these parameters do not affect the 
intervals of the other oscillation parameters. 


When combining the results for the T2K joint oscillation analyses in Secs. IX and |X| with 
the results from the reactor experiments, we use the weighted average of the results from the 
three reactor experiments Daya Bay, RENO, and Double Chooz which is: (sin 2 20i3) reactor = 
0.095±0.01 [lOOj. In terms of the parametrization that we use in this paper, (sin 2 di 3 ) reactor . = 
0.0243 ± 0.0026. 
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TABLE XIX: Predicted number of CC candidates and v e CC candidates for an 
exposure of 6.57 xlO 20 POT with and without oscillations and with oscillations using the 
typical parameter values: sin 2 6 Vi = 0.306, A/ 7 % = 7.5 x 10 _ 5 eV 2 /c 4 , sin 2 6*23 = 0.5, 
Amj 2 = 2.4 x 10 - 3 eV 2 /c 4 , sin 2 0 13 = 0.0243, Sep = 0 and normal mass hierarchy. The total 
numbers are broken down into the intrinsic beam components (those without an arrow) 

and oscillated components. 



Vf, 

Osc. 

CC 

No osc. 

v e 

Osc. 

CC 

No osc. 

u. 

116.46 

431.77 

0.94 

1.38 

v e -» Vn 

0.16 

0 

0.00 

0 

iV 

7.81 

13.92 

0.05 

0.06 

V e 

0.26 

0.27 

3.13 

3.38 

V e 

0.26 

0 

16.55 

0 

Ve 

0.02 

0.02 

0.15 

0.16 

Vfj, V e 

0.00 

0 

0.22 

0 

Total 

124.98 

445.98 

21.06 

4.97 


C. Oscillation parameter estimation 


Sections |VIII| - |X| describe analyses which use T2K and external data to estimate oscillation 
parameters and provide frequentist confidence intervals or Bayesian credible intervals. Using 
the disappearance channel alone, the atmospheric oscillation parameters are studied using 
frequentist approaches. The disappearance and appearance channels are used in combination 
to study a larger set of oscillation parameters, using frequentist and Bayesian approaches. 
This section describes general methods that are applied in these analyses. 

The oscillation analyses compare the event rate and distribution of the reconstructed 
neutrino energies for the observed v^ CC and v e CC candidate events recorded by the far 
detector, selected as described in Sec. |VI A[ with model predictions. The overall number of 


predicted events for typical oscillation parameter values and without oscillations are shown 
in Tab. ITTXl 

Point estimates for the oscillation parameters are those that maximize a likelihood func- 
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tion (or the posterior probability density for Bayesian analyses) that accounts for T2K-SK 
data, as well as internal control samples and external data. The observed numbers of events 
in SK are treated as outcomes of Poisson distributions. Systematic uncertainties are encap¬ 
sulated by the systematic parameters and their covariance matrices, defined in Sec.s 0Y1 


These provide a convenient mechanism to connect the separate analyses of the neutrino 
beamline, neutrino interactions, near detectors, and far detector to the full oscillation anal¬ 
yses. The analyses use different approaches to deal with the large number of oscillation and 
nuisance parameters, and report intervals based on either frequentist or Bayesian methods. 

With the large number of oscillation and nuisance parameters involved, it is not possible 
to calculate confidence intervals for a subset of the parameters with a method that guaran¬ 
tees frequentist coverag^j] for any possible values of the remaining parameters. Instead, a 
pragmatic approach is followed by reducing the high dimensionality of the likelihood func¬ 
tions through either profiling or marginalization. The profile likelihood, a function of only 
the subset of parameters of interest, is the likelihood maximized over the remaining param¬ 
eters. The marginal likelihood is found by integrating the product of the likelihood function 
and priors over all parameters, except those of interest. In the case of linear parameter 
dependence and where the nuisance parameters appear in a Gaussian form, the profile and 
marginal likelihood functions will be identical and can be used to produce intervals with 
correct frequentist coverage. For the neutrino oscillation analysis, the parameter depen¬ 
dence is non-linear and as a result the profile and marginal likelihoods differ, and frequentist 
coverage is not guaranteed. 

When practical, we use the Neyman approach of constructing a% confidence intervals 
whereby, for any value of the parameter(s) of interest, a% of possible data outcomes are 
accepted on the basis of a statistic. In our analyses, they are accepted if the likelihood 
ratio is larger than a critical value. The confidence interval is the set of all values for the 
parameter(s) for which the data are accepted. When physical boundaries or non-linearities 
appear in the parametrization, as in the case for the oscillation parameters, they can cause 
confidence intervals to be empty or misleadingly small. In order to reduce the chance of 
producing such confidence intervals, we use the likelihood ratio recommended by Feldman 
and Cousins [ 101J to form the interval. When producing joint intervals for two oscillation 
1 coverage demands that in an ensemble of repeated experiments, a% of the a% confidence intervals contain 
the true parameter(s). Coverage in presence of systematic uncertainty is difficult to define, in part due to 


the definition of an appropriate ensemble. 
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parameters, this approach is not always computationally practical, and instead approximate 
intervals are shown using contours of the likelihood ratio, sometimes referred to as the 
constant Ay 2 method. 


To construct Bayesian credible intervals, the posterior probability density function of the 
oscillation and nuisance parameters is calculated as the product of the likelihood function for 
the SK data with prior probability functions for the parameters. The Markov-Chain Monte 
Carlo (MCMC) method using the Metropolis-Hastings algorithm [102] is used to efficiently 
produce a set of points that populate the full parameter space proportional to the posterior 
probability density function. The chain is the set of accepted stepping points in a random 
walk through parameter space, in which a proposed step from point A to a point B with 
lower density is accepted with a probability equal to the ratio of the densities /(£>)/ f(A) and 
is always accepted when the density increases. When a step is not accepted, the last point 
in the chain is repeated, and another random step from that point is proposed. With the 
chain, consisting typically of millions of points, a% highest-posterior-density (HPD) credible 
intervals umi are constructed by selecting the region of highest density that contain a% 
of all the points. HPD intervals are constructed such that no point in parameter space 
outside the interval has a higher probability density than any point inside the interval. This 
is done for one or two parameters of interest, and the values of the remaining parameters 
are ignored in the process, equivalent to producing a set of points distributed according to 
the marginalized posterior probability density function. Unlike the frequentist approaches 
used, for which coverage is approximate, there are no approximations necessary to produce 
the credible intervals. 


The prior probability densities are, by default, uniform for the oscillation parameters over 
a large bounded region in the standard oscillation parametrization (|Am 2 |, sin 2 # 2 3 , sin 2 0 13 , 
Sqp ), multidimensional Gaussians for the nuisance parameters, and the prior probabilities 
for the two mass hierarchies are set to 0.5. As a result, the posterior probability density is 
proportional to the likelihood functions used for the frequentist analyses. Checks are made 
for alternative priors which are uniform in the oscillation angles, and the resulting interval 
boundaries are not strongly affected. 
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D. Characterizing systematic uncertainty 


The systematic parameters considered for the oscillation analyses can be grouped into 
three different categories: i) SK flux parameters and cross section parameters in common 
with ND280, ii) independent cross section parameters and iii) SK efficiencies, final state and 
secondary interactions (FSI+SI) and photo-nuclear (PN) parameters. The first category 
includes the systematic uncertainties related to the neutrino flux at SK and some cross 
sections, which are constrained by the near detector data as explained in Sec. |V} The values 
and uncertainties of these parameters used in the oscillation analyses are summarized in 
Tabs. XII and XIII The independent cross section parameters, described in Sec. Ill, are 


related to the nuclear model, therefore independent between the near and far detector as 
they contain different nuclei, or those which are common between the near and far detector 
but for which the near detector is insensitive. Table IVTH in Sec. Mil summarizes the values 
and uncertainties of the independent cross section parameters used for the SK oscillation 
analyses. Finally, the far detector efficiencies and uncertainties on final state, secondary and 


photo-nuclear interactions are described in Sec. VI A covariance matrix is computed for 
the uncertainties in this group; however, the uncertainty on the SK reconstructed energy 
scale, estimated to be 2.4%, is not included in the calculation of the covariance matrix, but 
considered as an independent systematic parameter. 

The effects of the systematic uncertainties on the predicted event rate are summarized 


in Tab. XX for the typical values of the oscillation parameters. In this table, the effects are 
presented as percentage uncertainties computed by throwing 10 6 toy experiments, varying 
only the systematics in the selected category (fixing the rest to their nominal values) and 
finding the RMS/mean of the distribution of number of events. 

Figure [26] shows the total error envelope combining all systematic uncertainties, calculated 
as the RMS from 10 6 toy MC experiments generated with randomized systematic parameters, 
taking into account all correlations between them, with and without the constraint from the 
ND280 data, showing a clear reduction of the error envelope when the constraint is applied. 
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TABLE XX: Relative uncertainty (la) on the predicted rate of ry CC and v e CC 

candidate events. 


Source of uncertainty 

^ CC 

u e CC 

Flux and common cross sections 



(w/o ND280 constraint) 

21.7% 

26.0% 

(w ND280 constraint) 

2.7% 

3.2% 

Independent cross sections 

5.0% 

4.7% 

SK 

4.0% 

2.7% 

FSI+SI(+PN) 

3.0% 

2.5% 

Total 



(w/o ND280 constraint) 

23.5% 

26.8% 

(w ND280 constraint) 

7.7% 

6.8% 




FIG. 26: Total error envelopes for the reconstructed energy distributions of CC (left) 
and v e CC (right) candidate events, using typical oscillation parameter values, with and 

without the ND280 constraint applied. 
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VIII. I/ M -> Vp ANALYSIS 


T2K has published several measurements of muon neutrino disappearance vmm- 
These measurements were performed within the framework of the PMNS oscillation model 


described in Section |VII A| and provided best-fit estimates and frequentist confidence inter¬ 
vals for the values of the mixing parameter siii 2 0 2 3 and the mass-squared splitting Am| 2 
(Amy) in the case of the normal (inverted) hierarchy. Each successive measurement an¬ 
alyzed a larger dataset, and the most recent measurement provides the world’s strongest 
constraint on sin 2 d 23 nn This section gives a more detailed description of that analysis 
and the study of multi-nucleon effects. Reducing the uncertainty on the values of these two 
parameters is important for measuring CP violation in neutrino oscillations by T2K and 
other current and future experiments. Furthermore, precise measurements of sin 2 0 2 3 could 
constrain models of neutrino mass generation (lOTHl 12] . 


Method 


The ^-disappearance analysis is performed by comparing the rate and spectrum of recon¬ 
structed neutrino energies, Eq. 0 - in the CC candidate event sample with predictions 
calculated from Monte Carlo simulation. The predicted spectrum is calculated by applying 


the survival probability in Eq. (16) to a prediction for the unoscillated rate and spectrum. 


These predictions are derived from our models of the total expected neutrino flux at the 
detector (explained in Sec. 0 and the cross section predictions for neutrino-nucleus interac¬ 
tions on water (described in Sec. EJ , which are constrained by near detector data (described 
in Sec. 0- and a GEANT3 model of particle interactions and transport in the SK detec¬ 
tor. The models of the flux, interaction physics, and detector include systematic parameters, 
whose uncertainties are accounted for in the analysis by using their corresponding covariance 
matrices. 

The oscillation parameters are estimated using two independent maximum likelihood fits 
to the reconstructed energy spectrum. The fits use different likelihoods and software in order 
to serve as cross-checks to each other. One analysis uses an extended unbinned likelihood 
(Ml), while the other uses a binned likelihood (M2). The log-likelihood definitions, ignoring 
constant terms, are: 
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• Ml likelihood , 


N d 


-2 In C(G, g,x s ,s) = - 2 ^ In f(E r ™\9, g, x s ,. 


i= 1 


+ 2 (^/V p (6>, <?, f s , s) - iV d In N p (9, g, x a , s) 

I A,?Tt/-1a z! i A Tir-tr-l 


( 19 ) 


A g 1 V- l Ag + Ax‘, VJAx, + As 1 V) 'A.?, 


and 


• M2 likelihood , 


N bi ns 

-2 In £( 0 , g,x s ,s) = - 2^2 N j ln 9,x s ,s) 

3 =1 

+ 27V P (0, g, x 9 , s) 

+ A/V"^ + AxjV^Ax, + As^I^As • 


( 20 ) 


In both definitions, 7V d and 1V P are the total number of data and predicted events re¬ 


spectively; 9 represents a vector of the PMNS oscillation parameters (Sec. VIIA); g is a 
vector containing the values of the systematic parameters constrained by the near detector 
(Tabs. XII|XIII ), x s are the cross section parameters not constrained by the near detector 
(Tab, |vh|), and s are the SK detector systematic parameters (Sec. |VI C[ ). A designates the 
difference between the systematic parameters and their nominal values, and V designates 
the covariance for the systematic parameters. 

For the Ml likelihood, /(E£® c |0, g, x s , s) is the probability density of observing an event 
with reconstructed energy, E™f, given values for the oscillation and systematic parameters. 
The value of f\E r Jf\9^ g, x s , s) is calculated with a linear interpolation between the bins of 
a histogram of the normalized energy spectrum. For the M2 likelihood, the number of data 
and predicted events in the j th reconstructed energy bin, Nj and Nj respectively, are used 
instead. 

Both z/^-disappearance fits consider a total of 48 parameters: 6 oscillation parameters, 16 
flux parameters, 20 neutrino interaction parameters and 6 parameters related to the response 
of SK. In order to find the best-fit values and confidence intervals for sin 2 023 and |Am 2 |, the 
profiled likelihood is maximized. Separate fits are performed for the different neutrino mass 
hierarchy assumptions. 
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B. Determining Confidence Intervals 


As explained in Sec. VIIC the Neyman method with the approach recommended by 
Feldman and Cousins (FC) was used to calculate confidence intervals for the two oscillation 
parameters, sin 2 0 23 and A'm 2 2 (Am 2 3 ), for the normal (inverted) hierarchy. The constant- 
Ay 2 method does not provide correct coverage due to the physical boundary near sin 2 20 23 = 1 
and because of the non-linear parametrization. Critical values of the FC statistic were deter¬ 
mined on a fine grid of the two oscillation parameters of interest using 10,000 toy datasets at 
each point. Each toy dataset had a set of values of the systematic parameters sampled from 
a multi-dimensional Gaussian having means at the nominal values, and covariances V. Each 
oscillation parameter, sin 2 0i 2 , Am 21 , and sin 2 0i 3 , is sampled from a Gaussian with mean 


and sigma values listed in Sec. |VII B| The values of 5cp are sampled uniformly between 
— 7 r and +7r. The systematic parameters and these additional oscillation parameters are 
removed from the likelihood function by profiling. In order to calculate an interval of just 
one oscillation parameter (sin 2 d 23 or Am 2 ), we determine the critical values by marginalizing 
over the second oscillation parameter. The marginalization assumes that the probability is 
proportional to the likelihood using T2K data. 


C. Results 


Both the Ml and M2 analyses find the point estimates sin 2 0 23 = 0.514 and A m\ 2 = 2.51 x 
10“ 3 eV 2 /c 4 when assuming the normal mass hierarchy and sin 2 6 l 23 = 0.511 and Am 2 3 = 


2.48 x 10 3 eV 2 /c 4 when assuming the inverted mass hierarchy. Table XXI summarizes these 


results from the Ml and M2 analyses. Likewise, the confidence intervals produced by Ml 
and M2 are similar. Since the Ml and M2 analyses are consistent with each other, only 


results from Ml are given below. Figure 27 shows the best-fit values of the oscillation 
parameters, the 2D confidence intervals calculated using the Feldman and Cousins method, 
assuming normal and inverted hierarchy, and the sensitivity at the current exposure. The 
size of the confidence interval found by the fit to the data is smaller than the sensitivity. 
This arises because the best-fit point is at the physical boundary corresponding to maximum 
disappearance probability. The amount by which the region is smaller is not unusual in an 
ensemble of toy MC experiments produced under the assumption of maximal disappearance. 
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FIG. 27: The 68% (dashed) and 90% (solid) CL intervals for the Ml ^-disappearance 
analysis assuming normal and inverted mass hierarchies. The 90% CL sensitivity contour 
for the normal hierarchy is overlaid for comparison. 


TABLE XXL Summary of the point estimates from the two independent 3-flavor muon 
neutrino disappearance oscillation frequentist analyses. 


Analysis 

MH 

A?n 2 2 or Am, 2 3 

(10" 3 eV 2 /c 4 ) 

sin 2 023 

TVfl R& 
exp 

Ml 

NH 

2.51 

0.514 

121.4 

Ml 

IH 

2.48 

0.511 

121.4 

M2 

NH 

2.51 

0.514 

121.5 

M2 

IH 

2.48 

0.511 

121.4 


The best-fit spectrum from the normal hierarchy fit compared to the observed spectrum is 


shown in Figure 28, showing as well the ratio of the number of observed events to the 
predicted number of events with sin 2 023 = 0. The observed oscillation dip is significant and 
well fit by simulation. The calculated ID Feldman and Cousins confidence intervals are 


given in Table XXII Figure 29 shows the -2Aln£ distributions for sin 2 023 and |Am 2 | from 
the data, along with the 90% CL critical values. 
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FIG. 28: Top: Reconstructed neutrino energy spectrum for data, best-fit prediction, and 
unoscillated prediction. Bottom: Ratio of oscillated to unoscillated events as a function of 
neutrino energy for the data and the best-fit spectrum. 

TABLE XXII: 68% and 90% confidence level intervals for the ^-disappearance analysis. 



MH 

68% CL 

90% CL 

sin 2 023 

NH 

[0.458, 0.568] 

[0.428, 0.598] 

sin 2 023 

IH 

[0.456, 0.566] 

[0.427, 0.596] 

Am 2 2 (10 -3 eV 2 /c 4 ) 

NH 

[2.41, 2.61] 

[2.34, 2.68] 

Amf 3 (10 -3 eV 2 /c 4 ) 

IH 

[2.38, 2.58] 

[2.31, 2.64] 


D. Multi-Nucleon Effects Study 

Recently, experimental (67[ T13ffll5j and theoretical [231 25 , 1116!il29j results have sug¬ 
gested that the charged-current neutrino-nucleus scattering cross section at T2K energies 
could contain a significant multi-nucleon component. Such processes are known to be impor- 
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FIG. 29: Profiled -2AlnL as a function of sin 2 $23 (left) and |Am 2 | (right), for the normal 
and inverted mass hierarchy assumptions. The 90% CL critical values are indicated by the 

lines with points. 


tant in describing electron-nucleus scattering (for a review, see [130] ). but have not yet been 
included in the model of neutrino-nucleus interactions in our muon neutrino disappearance 
analyses. If such multi-nucleon effects are important, their omission could introduce a bias 
in the oscillation analyses. Since low energy nucleons are not detected in SK, such events 
can be selected in the QE sample and assigned incorrect neutrino energies. 

A Monte Carlo study was performed in order to explore the sensitivity of the analysis to 
multi-nucleon effects. The nominal interaction model includes pion-less delta decay (PDD), 
which can be considered to be a multi-nucleon effect. As an alternative, we turn off PDD 
and use a model by Nieves [23j to simulate multi-nucleon interactions for neutrino energies 
below 1.5 GeV. Pairs of toy Monte Carlo experiments including both near and far detector 
data were generated, one with the nominal and one with the alternative model. Each dataset 
in a pair was produced by using the same distribution of interacting neutrinos, in order to 
reduce statistical fluctuations in the comparison. Each pair of experiments used a different 
distribution of interacting neutrinos and a different set of systematic parameters sampled 
from multivariate Gaussian distributions. The complete analysis with near and far detector 
data is performed, assuming the nominal model in all cases. In so doing, the study properly 
accounts for the reduction in sensitivity to mis-modcling neutrino interactions when using 
near detector data to constrain flux and cross section parameters. The differences in the 
point estimates for the oscillation parameters for the two samples in each pair are shown in 
Figure [30| The overall bias for both is negligible, compared to the precision obtained for the 
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FIG. 30: Difference in the point estimates of sin 2 023 (left) and | Am 2 1 (right) between pairs 
of toy MC datasets with and without including multi-nucleon effects. 

parameters. However, the additional variation in sin 2 023 is about 3%, comparable to the size 
of other systematic uncertainties. The bias was evaluated at sin 2 # 2 3 = 0.45 to avoid the the 
physical boundary at maximal disappearance which could reduce the size of the apparent 
bias. For the present exposure, the effect can be ignored, but future analyses will need to 
incorporate multi-nucleon effects in their model of neutrino-nucleus interactions. 
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IX. JOINT u, x DISAPPEARANCE AND v e APPEARANCE ANALYSIS USING 
A FREQUENTIST APPROACH. 


This section describes the joint 3-flavor oscillation analysis performed by combining the 
Vn disappearance and v e appearance channels using a frequentist approach. The oscillation 


parameters, 9 = |Am 2 |, sin 2 # 2 3 , shr#i 3 , and 5cp, described in Sec. VII A, are simultaneously 


determined. This is done by comparing the reconstructed energy spectra of the z/„ CC and 


u e CC candidate events observed at SK, selected as described in Sec. [VlJ with the predicted 
reconstructed energy spectra. Point estimates of the oscillation parameters are found by 
minimizing the negative log-likelihood 

^\i bins 

X 2 = -2 In £(0, g, x s ,s) = 2N?(9, g, x s , s) - 2 ^ In N^(9, g, x s , s) 


i =1 

bins 


+ 2 Nf(6, g, x„ i)- 2 ^ JVj In JV* (0, g, x„ s) 


( 21 ) 


i =1 
zT-ir—l , 


A f V~ 1 Ag + A xi V~ s Ax s + As V~ L As . 


where (AC) is the observed number of CC ( u e CC) candidate events in the i th re¬ 
constructed energy bin, and A r V (NA) is the corresponding predicted number of events, 
calculated as a function of the oscillation parameters 6 and the vectors of systematic pa¬ 


rameters, g,x s ,s, as described for Eq. (20). 


The negative log-likelihood function is minimized using MINUIT. As explained in 


Sec. |VII[ the solar oscillation parameters are kept fixed for this analysis. To combine 
our measurement with the reactor measurements, we add the term, 

2 

. ( 22 ) 


2 i siiU 0 i 3 (sin 9\^j reac f 0r 

X reactor 


& reactor 


where (sin 2 6 * 13 ) reactor and (J reac tor are given in Sec. VIIB 


When maximizing the likelihood, the systematic parameters are allowed to vary in a 
wide range [-5er, +5cr] (where a is the square root of the corresponding diagonal element 
in the covariance matrix), with the exception of the spectral function parameter which 
is constrained to lie between 0 (RFG) and 1 (SF). A total of 64 systematic parameters, 
representing uncertainties in the far detector efficiencies, the reconstructed neutrino energy 
scale, final state and secondary interactions, the flux prediction, and the relevant neutrino 
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TABLE XXIII: Point estimates of the oscillation parameters for the joint 3-flavor 

oscillation frequentist analysis. 


MH 

Am | 2 or Ara 2 3 
(10 -3 eV 2 /c 4 ) 

sin 2 #23 

sin 2 #13 

5cp 

1X1 err p 

N 1 Re 

exp 

Ay 2 

NH 

2.51 

0.524 

0.0422 

1.91 

119.9 

28.00 

0.01 

IH 

2.49 

0.523 

0.0491 

1.01 

119.9 

28.00 

0.00 


interaction models, are considered. As with the disappearance analyses, the fit to the ND280 
near detector data described in Sec. [V] is applied as a multivariate Gaussian penalty term 
to constrain the flux uncertainties and cross sections common to the near and far detectors. 

The 1-dimensional limits and 2-dimensional confidence regions reported in this analysis 
are constructed using the constant Ay 2 method [99J with respect to a 4-dimensional best-fit 


point obtained by minimizing Eq. (21). An exception is the (sin 2 #i 3 , &cp) space without the 
reactor measurement, as that analysis has little power to constrain Sep- For that case, a 
best-fit value of sin 2 # 13 is found for fixed values of 5 C p in the interval [-7T, 7 t] (divided into 51 
bins), resulting in 1 -dimensional confidence regions for different values of 5qp with respect to 
a line of best-fit points. For the T2K data fit combined with the reactor constraint, described 


in Sec. IX B the Feldman and Cousins method [101] is used to produce confidence intervals, 
by finding critical values of Ay 2 as a function of 5cp and we report excluded regions for 

&CP- 


Results 


Point estimates for the oscillation parameters and the expected number of events are 
summarized in Tab. XXIII Notably, the value obtained for sin 2 #i 3 by T2K is larger than 


the value found by the reactor experiments, the best-fit value of sin 2 #23 is consistent with 
maximal disappearance, and the difference in Ay 2 between the solutions for each mass 
hierarchy is negligible. 

The profiled Ay 2 of each oscillation parameter was obtained by minimizing the negative 
log-likelihood with respect to the systematic parameters and other three oscillation param¬ 


eters using MINUIT. Figure 31 presents the profiled Ay 2 of each oscillation parameter, 
















FIG. 31: Profiled Ay 2 for the joint 3-flavor oscillation analysis without using reactor data. 
The parameter |Am 2 represents A m\ 2 or Afor normal and inverted mass hierarchy 
assumptions respectively. The horizontal lines show the critical Ay 2 values for one 
dimensional fits at the 68 % and 90 % CL (Ay 2 = 1.00 and 2.71 respectively). 

comparing the results for the normal and the inverted mass hierarchy. From these figures, 
the la intervals estimated using the Ay 2 = 1 criterion are: 

sin 2 6*23 = 0.5241°;°" (NH) sin 2 0 23 = 0.523^;°^ (IH) 
sin 2 033 = 0.0421K (NH) sin 2 d 13 = 0.0491°$* (IH) 

Am 2 2 = 2.51^0^2 (10~ 3 eV 2 /c 4 , NH) Am 2 3 = 2.49±g;i| (10- 3 eV 2 /c 4 , IH). 

Figure [32] presents the 68 % and 90% CL regions for the two mass hierarchy assumptions 
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FIG. 32: 68 % (dashed) and 90% (solid) CL regions, from the analysis without using 
reactor data, with different mass hierarchy assumptions using Ay 2 with respect to the 
best-fit point - that from the inverted hierarchy. The parameter |Am 2 | represents A m\ 2 or 
Am 2 3 for normal and inverted mass hierarchy assumptions respectively. The lower left plot 
shows ID confidence intervals in sin 2 $13 for different values of Scp- 


in the four 2 -dimensional oscillation parameter spaces (sin 2 $ 23 , Aml 2 ), (sin 2 0i 3 , Am? 3 ), 
(sin 2 #i 3 , Scp ), and (sin 2 0 2 3 , sin 2 #i 3 ), constructed using constant Ay 2 with respect to the 
inverted hierarchy best-fit point. 























TABLE XXIV: Point estimates of the oscillation parameters for the joint 3-flavor 
oscillation frequentist analysis combined with the results from reactor experiments. 


MH 

A m | 2 or Ara 2 3 
(10 -3 eV 2 /c 4 ) 

sin 2 #23 

sin 2 #13 

5cp 

TVfl Rf 1 

1X1 err p 

N 1 Re 

exp 

Ay 2 

NH 

2.51 

0.527 

0.0248 

-1.55 

120.4 

25.87 

0.00 

IH 

2.48 

0.533 

0.0252 

-1.56 

121.2 

23.57 

0.86 


B. Results for T2K combined with the reactor experiment result 


The point estimates for the oscillation parameters and the predicted number of events, 
when the reactor measurements are included in the likelihood function, are given in 


Tab. XXIV The estimate for sin # 13 is smaller than the result obtained with T2K data 


only, shown in Tab. XXIII The likelihood is maximum for normal mass hierarchy and for 


5cp = — 7t/ 2, where the appearance probability is largest, as shown in Fig. 25 

The profiled Ay 2 as a function of each oscillation parameter are presented in Fig. 33 
and the 68% and 90% CL regions for the two mass hierarchies constructed using Ay 2 with 


respect to the best-fit point, the one for the normal hierarchy, are presented in Figs. 34 and 


The confidence regions obtained in the (sin 2 # 2 3 , |Am 2 |) space are compared with the 


results from Super-Kamiokande [131] and the MINOS [132] experiments in Fig. 36 The 
results from T2K and MINOS used the latest value of sin 2 #i 3 from [ 100 ] to fit this parameter 
whereas the result from SK has sin 2 #i 3 fixed to the previous reactor value in |99j . In the 
three analyses 5cp was removed by profiling. 

An analysis using the Feldman and Cousins method was performed for the measurement 
of 5cp including a reactor constraint by creating 4000 toy MC experiments at fixed values of 
5cp in the interval [-7r, 7r] (divided into 51 bins), taking into account statistical fluctuations 
and systematic variations. The other three oscillation parameters are removed by profiling 
following the 3-dimensional Ay 2 surface obtained as a result of the joint fit with the reactor 
constraint. The values of the critical Ay 2 calculated using these toy experiments are overlaid 


with the curve of Ay 2 as a function of 5cp in Fig. 37, and give the following excluded regions 
for 5cp at the 90% C.L: 5cp= [0.15,0.83]7r for normal hierarchy and 5cp= [—0.08,1.09]7t for 
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FIG. 33: Profiled Ay 2 for the joint 3-flavor oscillation analysis combined with the results 
from reactor experiments. The parameter |Am 2 | represents A m\ 2 or Am 2 3 for normal and 
inverted mass hierarchy assumptions respectively. The horizontal lines show the critical 
Ay 2 values for one dimensional fits at the 68 % and 90 % CL (Ay 2 = 1.00 and 2.71 

respectively). 


inverted hierarchy. 

In order to thoroughly cross-check the analysis described above, an alternate frequentist 
joint fit analysis was performed which differs in the treatment of the systematic errors. This 
originated as part of an effort to simplify and reduce the computing power needed for the 
analysis and to perform a study of the future sensitivity of the experiment |133j . A new set 
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FIG. 34: 68% (dashed) and 90% (solid) CL regions from the analysis that includes results 
from reactor experiments with different mass hierarchy assumptions using Ay 2 with 
respect to the best-fit point, the one from the fit with normal hierarchy. The parameter 
|A/n 2 | represents Am^ or Am 2 3 for normal and inverted mass hierarchy assumptions 

respectively. 


of systematic parameters is used; they multiply the nominal expected number of or v e 
events, with one parameter for each reconstructed energy bin. Results from the alternate 


analysis agree with the results presented in Secs. IX A and IX B 
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FIG. 35: Comparison of 68% (dashed) and 90% (solid) CL regions combined with the 
results from reactor experiments with different mass hierarchy assumptions using Ay 2 with 
respect to the best-fit point, the one from the fit with normal hierarchy. The parameter 
| Am 2 1 represents Am| 2 or Am 2 3 for normal and inverted mass hierarchy assumptions 

respectively. 
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sin 2 9 73 


FIG. 36: 68% (dashed) and 90% (solid) CL regions for normal (top) and inverted (bottom) 
mass hierarchy combined with the results from reactor experiments in the (sin 2 # 23 , Am§ 2 ) 
space compared to the results from the Super-Kamiokande [131] and MINOS [132] 

experiments. 
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FIG. 37: Profiled A\ /2 as a function of Scp with the results of the critical A% 2 values for 
the normal and inverted hierarchies for the joint fit with reactor constraint, with the 

excluded regions found overlaid. 
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X. JOINT U tl -> Vp AND -> n e BAYESIAN ANALYSIS 


This section describes a complementary approach to the analysis detailed in Sec. IX 


which uses Bayesian techniques to extract most probable values of oscillation parameters 
and their uncertainties. Bayesian inference analysis methods construct posterior proba¬ 
bilities of a hypothesis given the data observed by combining prior information with the 
likelihood function. This technique allows one to naturally include prior information about 
systematic parameters and external experimental data in the interpretation of the results 
of the experiment. Another distinguishing feature for this analysis is the fact that full 
marginalization of systematic parameters is achieved intrinsically, without the assumption 
that the observables are linear functions of the systematic parameters, taking into account 
the actual dependencies on the nuisance parameters. 

The posterior distribution, produced using Bayes’ theorem, is too difficult to compute 
analytically. We use two numerical methods to perform the high-dimensional integral nec¬ 
essary when computing the posterior distribution: a Markov Chain Monte Carlo (MCMC) 
in Sec. |X A| and a sampling method in Sec. |X B| which is used as a cross-check. 


A. Joint Near-Far Markov Chain Monte Carlo Analysis 

1. Point estimates 

To extract information about the point estimate of oscillation parameters from the pos¬ 
terior distribution generated by the MCMC, the density of points in 4-dimensional space 
was estimated using a kernel density estimator (KDE) [13411135] , A KDE estimates a PDF 
by smearing the discrete points of a MCMC in the 4 dimensions of interest. The Gaussian 
width of the smearing was set to be variable, and inversely proportional to the local density 
of MCMC points; this technique counters potential under-smoothing in low density regions 
and potential over-smoothing in high density regions. The maximum of the PDF produced 
by the KDE was then maximized using MINUIT to find the most probable value. In the 
case of using only T2K data, there is little sensitivity to the Scp parameter, and so a line 
of most probable values was created by finding the 3-dimensional density of the MCMC at 
a series of values of 5 C p- 
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2. Samples 


Unlike the frequentist analyses described above, the joint near-far analysis does not use 
the covariance matrix produced by the ND280 analysis described in Sec. [Vj Instead, this 
analysis is performed simultaneously with the three ND280 CC samples, and the SK 
CC, and SK v e CC samples. By fitting all samples simultaneously, this analysis avoids 
any error coming from neglecting non-linear dependencies of the systematic parameters 
constrained by ND280 analysis on the oscillation parameters. 

The systematic uncertainties used for the ND280 samples are nearly identical to those in 
Sec. |V| with the following exceptions: the uncertainties on the cross section ratios c Ue /c 
and (Tp/cr^ are applied and the NC normalization uncertainties are divided into NCl 7 r°, 
NCItt 1 * 1 , NC coherent, and NCOther for all samples. Additionally, the number of bins in the 
ND280 detector systematic covariance matrix is reduced to 105, in order to reduce the total 
number of parameters. There are no differences in the systematic uncertainties for the SK 
samples. Ignoring constant terms, the negative log of the posterior probability is given by, 


N D280bins 

-In (P)= E ^(S,x,d)-N?tnN?(b,xJ) 

i 

Nn, bins 

i 


Ne bins 

i 

+ \Ai r V b 'Ab+\Ax T V-'Ax+\A3 r V d ;'Ad 
+ 1 Af'V-'AS+ \ L . 


(23) 


The vector 9 sr contains the solar oscillation parameters and for combined fits with reac¬ 
tor data sin 2 29 13 , with priors described in Sec. VIIB The priors on the other oscillation 
parameters of interest are uniform in sin 2 0 13 between 0 and 1 , sin 2 6*23 between 0 and 1 , 
|AmJ 2 | between 0.001 and 0.005 eV 2 /c 4 , and Sep between — n and 7 r. Additionally, the 
prior probability of the normal hierarchy and inverted hierarchy are each 0.5. Priors for the 
systematic parameters are the multivariate Gaussian terms shown, with the exception of the 
cross section spectral function parameters which are given a uniform prior between 0 and 1 . 

In this analysis, both ND280 and SK MC sample events are weighted individually for 
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TABLE XXV: Most probable values for oscillation parameters from Bayesian analysis. 



Hierarchy 

Am .32 

sin 2 $23 

sin 2 #13 

Scp 

Analysis 


10 ~ 3 eV 2 /c 4 




T2K-only 

Inverted 

2.571 

0.520 

0.0454 

0 (fixed) 

T2K+reactor 

Normal 

2.509 

0.528 

0.0250 

-1.601 


TABLE XXVI: 68 % Bayesian credible intervals for oscillation parameters. 



|Am| 2 | 

sin 2 #23 

sin 2 #13 

Analysis 

10 " 3 eV 2 /c 4 



T2K-only 

[2.46, 2.68] 

[0.470, 0.565] 

[0.0314 ,0.0664] 

T2K+reactor 

[2.40, 2.62] 

[0.490, 0.583] 

[0.0224, 0.0276] 


all parameters in the analysis. This means that each PDF is rebuilt from the MC at every 
iteration of the MCMC. This has the advantage of retaining shape information within each 
bin of the PDF, especially desirable for the oscillation parameters, and also allows a more 
natural treatment of certain parameters such as the SK energy scale uncertainty which 
may cause events to migrate between bins. The increase in computational load was offset 
by performing certain calculations on GPUs, including the event-by-event calculation of 
oscillation probability [136] . 


3. Results 


The MCMC was run with 5.6 x 10' steps using only T2K data, and for 1.4 x 10 8 steps 
for T2K data combined with reactor experiment results. The most probable values for the 


oscillation parameters for both analyses are shown in Table XXV For the T2K-only analysis, 


the values are shown for Scp= 0, as the analysis has little sensitivity to the value of Scp- The 
ID credible intervals, marginalized over all other parameters, including mass hierarchy, 


for each of the parameters except Scp are shown in Tabic XXVI 


Figures [38] and |39] show the Scp versus sin 2 #13 and A m\ 2 versus sin 2 #23 credible regions for 


the T2K-only and T2K+reactor analyses. Note that the contours in Fig. 38 are marginalized 
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- T2K+Reactor 68% Credible Region - T2K Only 68% Credible Region 

- T2K+Reactor 90% Credible Region - T2K Only 90% Credible Region 

• T2K+Reactor Best Fit Point T2K Only Best Fit Line 


FIG. 38: Credible regions for sin 2 #13 and Scp for T2K-only and T2K+reactor combined 
analyses. These are constructed by marginalizing over both mass hierarchies. For the 
T2K-only analysis, the best fit line is shown instead of the best fit point because the 

analysis has little sensitivity to Scp- 


over the mass hierarchy; in particular, the most probable value line appears to be offset 
from the center of the credible region. This is because the most probable value line is for 
the preferred inverted hierarchy, and the credible intervals are marginalized over hierarchy. 


Fig. 40 shows the posterior probability for Scp with 68 % and 90% credible intervals for the 


T2K+reactor combined analysis. Figure 41 shows comparisons of SK CC and v e CC 
candidate events with the best-fit spectra produced from the T2K-only and T2K+reactor 
combined analyses. Each best-fit spectrum is formed by calculating the most probable value 
for the predicted number of events in each energy bin, using all of the MCMC points from 
the corresponding analysis. The fit spectrum for CC events does not change appreciably 
when the reactor prior is included, but the v e CC fit spectrum shows a noticeable reduction 
in the number of events. 
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FIG. 39: Credible regions for sin 2 $23 and Am | 2 for T2K-only and T2K+reactor combined 
analyses. The normal hierarchy corresponds to positive values of A m\ 2 and the inverted 

hierarchy to negative values. 


Figures [42] and |43] show the posterior PDFs for the oscillation parameters both singly 
and pairwise, using MCMC points from the inverted and normal hierarchy respectively, 
which reflect the most probable mass hierarchy for the T2K-only and T2K+reactor analysis 
respectively. The plots along the diagonal show the posterior PDFs for each of the four 
oscillation parameters of interest, marginalized over all other parameters, except for the 
mass hierarchy. The off-diagonal elements show the pairwise posterior PDFs. 

Another interesting feature of this analysis is that it provides a natural way to study 
the preference of the data for normal versus inverted hierarchy and lower versus upper 
octant in 6 * 23 . This is done simply by comparing the total probability (that is, the number 


of MCMC steps) in the region of interest. Table XXVII shows the probability for the 
various cases for the T2K-only analysis. Note that the inverted hierarchy is preferred in 
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FIG. 40: The posterior probability for Scp, marginalized over all other parameters, 
including mass hierarchy, for the T2K+reactor combined analysis. 



FIG. 41: T2K-only and T2K+reactor prior best-fit spectra overlaid with SK CC and v e 

CC candidate samples. 


this analysis, but the posterior odds raticj^] is only 1.2. Table XXVIII shows the same for 
the T2K+reactor combined analysis. In this analysis, the normal hierarchy is preferred, but 
with a posterior odds ratio of 2.2, the inverted hierarchy is not significantly excluded with 
the present analysis. 


with the prior odds assumed to be 1, the posterior odds ratio is equivalent to the Bayes Factor 
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FIG. 42: Distributions of posterior probability between the oscillation parameters of 
interest for the T2K-only analysis. These posteriors use only MCMC points that are in the 

inverted hierarchy. 


TABLE XXVII: Model comparison probabilities for normal and inverted mass hierarchies, 
as well as upper and lower octants, without including reactor data. 



NH 

IH 

Sum 

sin 2 023 < 0.5 

0.165 

0.200 

0.365 

sin 2 023 > 0.5 

0.288 

0.347 

0.635 

Sum 

0.453 

0.547 

1.0 


To evaluate the dependency of this analysis on the form of the prior of the oscillation 
parameters, the analysis was repeated with a uniform prior in 9 13 and 623 . The credible in¬ 
tervals and model comparison probabilities do not change appreciably with these alternative 
priors. 
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FIG. 43: Distributions of posterior probability between the oscillation parameters of 
interest for the T2K+reactor analysis. These posteriors use only MCMC points that are in 


the normal hierarchy. In comparing to Fig. 42, note the change in scales for some 

parameters. 


B. Cross-check analysis 

A second Bayesian joint analysis (JB2) is used to cross-check the results from the analysis 
described above (JB1). Like the frequentist analyses, JB2 uses the output from the ND280 
analysis described in Sec. [V]to constrain some of the systematic uncertainties, by applying 
them as prior probability densities. Also, JB2 does not use (by default) the reconstructed 
energy spectrum for v e candidate events, but instead the 2D distribution of the momentum 
and angle with respect to beam direction (p e , 6 e ) of the particle reconstructed as an electron 
in those events. This is similar to what was used in the previously reported electron neutrino 
appearance observation (4j. JB2 can also use the shape of the reconstructed energy spectrum 
for v e candidate events, so that the results of the two analyses can be compared in both 
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TABLE XXVIII: Model comparison probabilities for normal and inverted mass hierarchies, 
as well as upper and lower octants, including reactor data. 



NH 

IH 

Sum 

sin 2 #23 < 0.5 

0.179 

0.078 

0.257 

sin 2 #23 > 0.5 

0.505 

0.238 

0.743 

Sum 

0.684 

0.316 

1.0 


cases. On a technical level, MCMC is not used in this second analysis to marginalize 
over the nuisance parameters; the integration is done numerically by averaging the posterior 
probability over 10,000 throws of those parameters following their prior distribution. Finally, 
a second technical difference is that in JB2 the weighting is not done event by event but by 
(Pe, 0 e ) bin. 


C. Comparison of analyses 


1. Comparison of Bayesian joint analyses 


The results obtained with the two joint Bayesian analyses are very similar, both in terms 
of posterior probabilities for the different models and credible intervals for the oscillation 
parameters. The comparison in the case of the posterior probability for Scp is shown in 


Fig. [44} the posterior probabilities obtained by the two analyses are similar, and most of the 
difference comes from JB2 using the (p e , 6 e ) spectrum shape for u e candidate events instead 
of the reconstructed energy spectrum shape as JB1 does. This also shows that at the current 
statistics, fitting the near and far detector samples at the same time and using the output 
of the near detector analysis described in Sec. [V]are equivalent. 


2. Treatment of the systematic uncertainties 


We also compare, using JB2, the marginalization and profiling approaches described in 
Sec. VTTp to reduce the dimensionality of the likelihood. In the case of Scp , the marginal 
(obtained by integrating the product of the likelihood and priors over the nuisance parame¬ 
ters) and profile (obtained by maximizing the likelihood with respect to those parameters) 
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FIG. 44: Posterior probabilities for §cp obtained by the two joint Bayesian analyses using 

the reactor experiments prior for sin 2 0 13 . 



FIG. 45: Marginal and profile likelihoods of the T2K data with reactor constraint 

assuming normal hierarchy. 


likelihoods are visibly different, as can be seen on figure 45 Such differences are expected 
as some of the nuisance parameters appear in a non-Gaussian form and have a non-linear 
dependence. Within the Bayesian framework, only marginalization is well motivated. 
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XI. CONCLUSIONS 


With the data collected between 2010 and 2013 we have analyzed the z/^-disappearance 
to estimate the two oscillation parameters, |Am 2 | and sin 2 0 2 3 - For the first time, we have 
used a combined analysis of ^-disappearance and z/ e -appearance, to advance our knowledge 
of the oscillation parameters |Ara 2 |, sin 2 # 23 , sin 2 0i3, 5cp, and the mass hierarchy. 

Uncertainty arising from systematic factors has been carefully assessed in the analyses and 
its effect is small compared to statistical errors. Our understanding of neutrino oscillation 
will continue to improve as we collect more data in the coming years, in both neutrino and 
anti-neutrino mode [133] . The general approach followed in this paper that couples the 
separate analysis of the beamline, neutrino interactions, near detectors, and far detector, 
through sets of systematic parameters and their covariances, will be extended to deal with 
additional information from anti-neutrino data and from additional selections with the near 
detector data. 
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